Index: /issm/trunk-jpl/src/m/archive/arch.py =================================================================== --- /issm/trunk-jpl/src/m/archive/arch.py (revision 21302) +++ /issm/trunk-jpl/src/m/archive/arch.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import math import struct @@ -150,16 +150,16 @@ def write_vector(fid,data): # {{{ """ - Procedure to write a numpy array to an arch file + Procedure to write a np.array to an arch file """ # Make sure our vector is the correct shape. # Reshape it into a row vector if it is not correct. if isinstance(data,(bool,int,long,float)): - data=numpy.array([data]) + data=np.array([data]) elif isinstance(data,(list,tuple)): - data=numpy.array(data).reshape(-1,1) - - if numpy.ndim(data) == 1: - if numpy.size(data): - data=data.reshape(numpy.size(data),1) + data=np.array(data).reshape(-1,) + + if np.ndim(data) == 1: + if np.size(data): + data=data.reshape(np.size(data),) else: data=data.reshape(0,0) @@ -216,5 +216,5 @@ rows=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0] cols=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0] - raw_data=numpy.zeros(shape=(rows,cols),dtype=float) + raw_data=np.zeros(shape=(rows,cols),dtype=float) for i in xrange(rows): raw_data[i,:]=struct.unpack('>%dd' % cols,fid.read(cols*struct.calcsize('>d'))) @@ -259,5 +259,5 @@ elif format.shape[0] == 1 and format.shape[1] == 1: code=2 - elif isinstance(format,(list,tuple,numpy.ndarray)): + elif isinstance(format,(list,tuple,np.ndarray)): code=3 else: Index: /issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py (revision 21302) +++ /issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py (revision 21303) @@ -1,4 +1,4 @@ import os -import numpy as np +import numpy as np def PattynSMB(md,Tf): """ Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 21302) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 21303) @@ -1,4 +1,4 @@ import os -import numpy +import numpy as np from ContourToMesh import ContourToMesh @@ -14,16 +14,16 @@ #node on Dirichlet - pos=numpy.nonzero(md.mesh.vertexonboundary) - md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices) + pos=np.nonzero(md.mesh.vertexonboundary) + md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) md.stressbalance.spcvx[pos]=0 md.stressbalance.spcvy[pos]=0 md.stressbalance.spcvz[pos]=0 - md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) - md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3)) + md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) + md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) #Dirichlet Values - if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: + if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: print " boundary conditions for stressbalance model: spc set as observed velocities" md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] @@ -39,18 +39,18 @@ #Deal with other boundary conditions - if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): - md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(md.balancethickness.thickening_rate)): + md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) print " no balancethickness.thickening_rate specified: values set as zero" - md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) - if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: - md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: + md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) if hasattr(md.mesh,'vertexonsurface'): - pos=numpy.nonzero(md.mesh.vertexonsurface)[0] + pos=np.nonzero(md.mesh.vertexonsurface)[0] md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface - if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices: - md.basalforcings.geothermalflux=50.*10**-3*numpy.ones((md.mesh.numberofvertices,1)) #50 mW/m^2 + if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices: + md.basalforcings.geothermalflux=50.*10**-3*np.ones((md.mesh.numberofvertices)) #50 mW/m^2 else: print " no thermal boundary conditions created: no observed temperature found" Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 21302) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 21303) @@ -1,4 +1,4 @@ import os -import numpy +import numpy as np from ContourToMesh import ContourToMesh import MatlabFuncs as m @@ -27,18 +27,18 @@ raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) - nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) + nodeonicefront=np.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) else: - nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool) + nodeonicefront=np.zeros((md.mesh.numberofvertices),bool) # pos=find(md.mesh.vertexonboundary & ~nodeonicefront); - pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0] - md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) - md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3)) + pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(nodeonicefront)))[0] + md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) + md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) #Icefront position - pos=numpy.nonzero(nodeonicefront)[0] + pos=np.nonzero(nodeonicefront)[0] md.mask.ice_levelset[pos]=0 @@ -53,10 +53,10 @@ values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1] segmentsfront=1-values - numpy.sum(segmentsfront,axis=1)!=numbernodesfront - segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0] + np.sum(segmentsfront,axis=1)!=numbernodesfront + segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0] #Find all nodes for these segments and spc them pos=md.mesh.segments[segments,0:-1]-1 else: - pos=numpy.nonzero(md.mesh.vertexonboundary)[0] + pos=np.nonzero(md.mesh.vertexonboundary)[0] md.stressbalance.spcvx[pos]=0 md.stressbalance.spcvy[pos]=0 @@ -64,10 +64,10 @@ #Dirichlet Values - if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: + if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: #reshape to rank-2 if necessary to match spc arrays - if numpy.ndim(md.inversion.vx_obs)==1: - md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,1) - if numpy.ndim(md.inversion.vy_obs)==1: - md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,1) + if np.ndim(md.inversion.vx_obs)==1: + md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,) + if np.ndim(md.inversion.vy_obs)==1: + md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,) print " boundary conditions for stressbalance model: spc set as observed velocities" md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] @@ -81,18 +81,18 @@ #Deal with other boundary conditions - if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): - md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(md.balancethickness.thickening_rate)): + md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) print " no balancethickness.thickening_rate specified: values set as zero" - md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) - if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: - md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: + md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) if hasattr(md.mesh,'vertexonsurface'): - pos=numpy.nonzero(md.mesh.vertexonsurface)[0] + pos=np.nonzero(md.mesh.vertexonsurface)[0] md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface - if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: - md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) + if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: + md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices)) else: print " no thermal boundary conditions created: no observed temperature found" Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 21302) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 21303) @@ -1,4 +1,4 @@ import os -import numpy +import numpy as np from ContourToMesh import ContourToMesh import MatlabFuncs as m @@ -29,25 +29,25 @@ raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile) incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) - vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1)) + vertexonicefront=np.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1)) else: #Guess where the ice front is - vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,1)) - pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0] + vertexonfloatingice=np.zeros((md.mesh.numberofvertices)) + pos=np.nonzero(np.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0] vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1. - vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.) + vertexonicefront=np.logical_and(np.reshape(md.mesh.vertexonboundary,(-1,)),vertexonfloatingice>0.) # pos=find(md.mesh.vertexonboundary & ~vertexonicefront); - pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0] - if not numpy.size(pos): + pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(vertexonicefront)))[0] + if not np.size(pos): print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually." - md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices) - md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) - md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3)) + md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) + md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) + md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) #Position of ice front - pos=numpy.nonzero(vertexonicefront)[0] + pos=np.nonzero(vertexonicefront)[0] md.mask.ice_levelset[pos]=0 @@ -62,10 +62,10 @@ values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1] segmentsfront=1-values - numpy.sum(segmentsfront,axis=1)!=numbernodesfront - segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0] + np.sum(segmentsfront,axis=1)!=numbernodesfront + segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0] #Find all nodes for these segments and spc them pos=md.mesh.segments[segments,0:-1]-1 else: - pos=numpy.nonzero(md.mesh.vertexonboundary)[0] + pos=np.nonzero(md.mesh.vertexonboundary)[0] md.stressbalance.spcvx[pos]=0 md.stressbalance.spcvy[pos]=0 @@ -73,5 +73,5 @@ #Dirichlet Values - if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: + if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: print " boundary conditions for stressbalance model: spc set as observed velocities" md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] @@ -80,6 +80,6 @@ print " boundary conditions for stressbalance model: spc set as zero" - md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2)) - pos=numpy.nonzero(md.mesh.vertexonboundary)[0] + md.hydrology.spcwatercolumn=np.zeros((md.mesh.numberofvertices,2)) + pos=np.nonzero(md.mesh.vertexonboundary)[0] md.hydrology.spcwatercolumn[pos,0]=1 @@ -89,20 +89,20 @@ #Deal with other boundary conditions - if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): - md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(md.balancethickness.thickening_rate)): + md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) print " no balancethickness.thickening_rate specified: values set as zero" - md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) + md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) - if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: - md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: + md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) if hasattr(md.mesh,'vertexonsurface'): - pos=numpy.nonzero(md.mesh.vertexonsurface)[0] + pos=np.nonzero(md.mesh.vertexonsurface)[0] md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface - if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: - md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) - md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3 #50mW/m2 + if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: + md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices)) + md.basalforcings.geothermalflux[np.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3 #50mW/m2 else: print " no thermal boundary conditions created: no observed temperature found" Index: /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py (revision 21302) +++ /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py (revision 21303) @@ -1,5 +1,5 @@ from MatlabFuncs import * from model import * -from numpy import * +from np.import * def love_numbers(value,*varargin): Index: /issm/trunk-jpl/src/m/classes/SMBcomponents.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 21303) @@ -38,14 +38,14 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.accumulation)): - self.accumulation=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.accumulation)): + self.accumulation=np.zeros((md.mesh.numberofvertices)) print " no SMB.accumulation specified: values set as zero" - if numpy.all(numpy.isnan(self.runoff)): - self.runoff=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.runoff)): + self.runoff=np.zeros((md.mesh.numberofvertices)) print " no SMB.runoff specified: values set as zero" - if numpy.all(numpy.isnan(self.evaporation)): - self.evaporation=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.evaporation)): + self.evaporation=np.zeros((md.mesh.numberofvertices)) print " no SMB.evaporation specified: values set as zero" Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -64,10 +64,10 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.s0p)): - self.s0p=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.s0p)): + self.s0p=np.zeros((md.mesh.numberofvertices)) print " no SMBd18opdd.s0p specified: values set as zero" - if numpy.all(numpy.isnan(self.s0t)): - self.s0t=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.s0t)): + self.s0t=np.zeros((md.mesh.numberofvertices)) print " no SMBd18opdd.s0t specified: values set as zero" @@ -90,6 +90,6 @@ if 'MasstransportAnalysis' in analyses: md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',[1]) - md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) - md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',[1]) md = checkfield(md,'fieldname','smb.rlapslgm','>=',0,'numel',[1]) @@ -98,5 +98,5 @@ md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',[1]) Index: /issm/trunk-jpl/src/m/classes/SMBforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBforcing.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/SMBforcing.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -33,6 +33,6 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.mass_balance)): - self.mass_balance=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.mass_balance)): + self.mass_balance=np.zeros((md.mesh.numberofvertices)) print " no SMBforcing.mass_balance specified: values set as zero" Index: /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py (revision 21303) @@ -40,18 +40,18 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.accumulation)): - self.accumulation=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.accumulation)): + self.accumulation=np.zeros((md.mesh.numberofvertices)) print " no SMB.accumulation specified: values set as zero" - if numpy.all(numpy.isnan(self.evaporation)): - self.evaporation=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.evaporation)): + self.evaporation=np.zeros((md.mesh.numberofvertices)) print " no SMB.evaporation specified: values set as zero" - if numpy.all(numpy.isnan(self.melt)): - self.melt=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.melt)): + self.melt=np.zeros((md.mesh.numberofvertices)) print " no SMB.melt specified: values set as zero" - if numpy.all(numpy.isnan(self.refreeze)): - self.refreeze=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.refreeze)): + self.refreeze=np.zeros((md.mesh.numberofvertices)) print " no SMB.refreeze specified: values set as zero" Index: /issm/trunk-jpl/src/m/classes/SMBpdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -94,10 +94,10 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.s0p)): - self.s0p=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.s0p)): + self.s0p=np.zeros((md.mesh.numberofvertices)) print " no SMBpdd.s0p specified: values set as zero" - if numpy.all(numpy.isnan(self.s0t)): - self.s0t=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.s0t)): + self.s0t=np.zeros((md.mesh.numberofvertices)) print " no SMBpdd.s0t specified: values set as zero" @@ -119,6 +119,6 @@ if 'MasstransportAnalysis' in analyses: md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',[1]) - md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) - md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',[1]) md = checkfield(md,'fieldname','smb.rlapslgm','>=',0,'numel',[1]) @@ -128,12 +128,12 @@ md = checkfield(md,'fieldname','smb.precipitation','NaN',1,'Inf',1,'timeseries',1) elif self.isdelta18o: - md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) - md = checkfield(md,'fieldname','smb.delta18o_surface','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.delta18o_surface','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) - md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) elif self.ismungsm: md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) @@ -141,7 +141,7 @@ md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','smb.Pfac','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) - md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) - md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,numpy.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.Pfac','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) + md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1) Index: /issm/trunk-jpl/src/m/classes/adinversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/adinversion.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/adinversion.py (revision 21303) @@ -12,15 +12,10 @@ from MatlabFuncs import * -from numpy import * +from np.import * # ADINVERSION class definition - # - # Usage: - # adinversion=adinversion(); - - class adinversion: @@ -75,6 +70,6 @@ - num_controls=numpy.numel(md.inversion.control_parameters) - num_costfunc=numpy.size(md.inversion.cost_functions,2) + num_controls=np.numel(md.inversion.control_parameters) + num_costfunc=np.size(md.inversion.cost_functions,2) @@ -95,12 +90,12 @@ if solution=='BalancethicknessSolution': - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1) - md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices, 1], float('Nan'),1) + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],float('Nan'),1) + md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices], float('Nan'),1) elif solution=='BalancethicknessSoftSolution': - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1) + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],float('Nan'),1) else: - md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1) - if not numpy.strcmp(domaintype(md.mesh),'2Dvertical'): - md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1) + md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],float('Nan'),1) + if not np.strcmp(domaintype(md.mesh),'2Dvertical'): + md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],float('Nan'),1) return md @@ -164,34 +159,34 @@ #process control parameters - num_control_parameters = numpy.numel(self.control_parameters); + num_control_parameters = np.numel(self.control_parameters); WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray'); WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer'); #process cost functions - num_cost_functions=numpy.size(self.cost_functions,2); + num_cost_functions=np.size(self.cost_functions,2); data=copy.deepcopy(self.cost_functions) - data[numpy.nonzero(self.cost_functions==101)] =['SurfaceAbsVelMisfit']; - data[numpy.nonzero(self.cost_functions==102)]=['SurfaceRelVelMisfit']; - data[numpy.nonzero(self.cost_functions==103)]=['SurfaceLogVelMisfit']; - data[numpy.nonzero(self.cost_functions==104)]=['SurfaceLogVxVyMisfit']; - data[numpy.nonzero(self.cost_functions==105)]=['SurfaceAverageVelMisfit']; - data[numpy.nonzero(self.cost_functions==201)]=['ThicknessAbsMisfit']; - data[numpy.nonzero(self.cost_functions==501)]=['DragCoefficientAbsGradient']; - data[numpy.nonzero(self.cost_functions==502)]=['RheologyBbarAbsGradient']; - data[numpy.nonzero(self.cost_functions==503)]=['ThicknessAbsGradient']; - data[numpy.nonzero(self.cost_functions==504)]=['ThicknessAlongGradient']; - data[numpy.nonzero(self.cost_functions==505)]=['ThicknessAcrossGradient']; - data[numpy.nonzero(self.cost_functions==506)]=['BalancethicknessMisfit']; - data[numpy.nonzero(self.cost_functions==601)]=['SurfaceAbsMisfit']; - data[numpy.nonzero(self.cost_functions==1001)]=['Outputdefinition1']; - data[numpy.nonzero(self.cost_functions==1002)]=['Outputdefinition2']; - data[numpy.nonzero(self.cost_functions==1003)]=['Outputdefinition3']; - data[numpy.nonzero(self.cost_functions==1004)]=['Outputdefinition4']; - data[numpy.nonzero(self.cost_functions==1005)]=['Outputdefinition5']; - data[numpy.nonzero(self.cost_functions==1006)]=['Outputdefinition6']; - data[numpy.nonzero(self.cost_functions==1007)]=['Outputdefinition7']; - data[numpy.nonzero(self.cost_functions==1008)]=['Outputdefinition8']; - data[numpy.nonzero(self.cost_functions==1009)]=['Outputdefinition8']; - data[numpy.nonzero(self.cost_functions==1010)]=['Outputdefinition10']; + data[np.nonzero(self.cost_functions==101)] =['SurfaceAbsVelMisfit']; + data[np.nonzero(self.cost_functions==102)]=['SurfaceRelVelMisfit']; + data[np.nonzero(self.cost_functions==103)]=['SurfaceLogVelMisfit']; + data[np.nonzero(self.cost_functions==104)]=['SurfaceLogVxVyMisfit']; + data[np.nonzero(self.cost_functions==105)]=['SurfaceAverageVelMisfit']; + data[np.nonzero(self.cost_functions==201)]=['ThicknessAbsMisfit']; + data[np.nonzero(self.cost_functions==501)]=['DragCoefficientAbsGradient']; + data[np.nonzero(self.cost_functions==502)]=['RheologyBbarAbsGradient']; + data[np.nonzero(self.cost_functions==503)]=['ThicknessAbsGradient']; + data[np.nonzero(self.cost_functions==504)]=['ThicknessAlongGradient']; + data[np.nonzero(self.cost_functions==505)]=['ThicknessAcrossGradient']; + data[np.nonzero(self.cost_functions==506)]=['BalancethicknessMisfit']; + data[np.nonzero(self.cost_functions==601)]=['SurfaceAbsMisfit']; + data[np.nonzero(self.cost_functions==1001)]=['Outputdefinition1']; + data[np.nonzero(self.cost_functions==1002)]=['Outputdefinition2']; + data[np.nonzero(self.cost_functions==1003)]=['Outputdefinition3']; + data[np.nonzero(self.cost_functions==1004)]=['Outputdefinition4']; + data[np.nonzero(self.cost_functions==1005)]=['Outputdefinition5']; + data[np.nonzero(self.cost_functions==1006)]=['Outputdefinition6']; + data[np.nonzero(self.cost_functions==1007)]=['Outputdefinition7']; + data[np.nonzero(self.cost_functions==1008)]=['Outputdefinition8']; + data[np.nonzero(self.cost_functions==1009)]=['Outputdefinition8']; + data[np.nonzero(self.cost_functions==1010)]=['Outputdefinition10']; WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray'); WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer'); Index: /issm/trunk-jpl/src/m/classes/autodiff.py =================================================================== --- /issm/trunk-jpl/src/m/classes/autodiff.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/autodiff.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from dependent import dependent from independent import independent @@ -104,6 +104,6 @@ if num_dependent_objects: names=[] - types=numpy.zeros(num_dependent_objects) - indices=numpy.zeros(num_dependent_objects) + types=np.zeros(num_dependent_objects) + indices=np.zeros(num_dependent_objects) for i,dep in enumerate(self.dependents): @@ -122,5 +122,5 @@ if num_independent_objects: names=[None] * num_independent_objects - types=numpy.zeros(num_independent_objects) + types=np.zeros(num_independent_objects) for i,indep in enumerate(self.independents): @@ -136,5 +136,5 @@ for indep in self.independents: - if not numpy.isnan(indep.fos_forward_index): + if not np.isnan(indep.fos_forward_index): index+=indep.fos_forward_index break @@ -153,5 +153,5 @@ for dep in self.dependents: - if not numpy.isnan(dep.fos_reverse_index): + if not np.isnan(dep.fos_reverse_index): index+=dep.fos_reverse_index break Index: /issm/trunk-jpl/src/m/classes/bamggeom.py =================================================================== --- /issm/trunk-jpl/src/m/classes/bamggeom.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/bamggeom.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np class bamggeom(object): @@ -10,12 +10,12 @@ def __init__(self,*args): # {{{ - self.Vertices=numpy.empty((0,3)) - self.Edges=numpy.empty((0,3)) - self.TangentAtEdges=numpy.empty((0,4)) - self.Corners=numpy.empty((0,1)) - self.RequiredVertices=numpy.empty((0,1)) - self.RequiredEdges=numpy.empty((0,1)) - self.CrackedEdges=numpy.empty((0,0)) - self.SubDomains=numpy.empty((0,4)) + self.Vertices=np.empty((0,3)) + self.Edges=np.empty((0,3)) + self.TangentAtEdges=np.empty((0,4)) + self.Corners=np.empty((0,1)) + self.RequiredVertices=np.empty((0,1)) + self.RequiredEdges=np.empty((0,1)) + self.CrackedEdges=np.empty((0,0)) + self.SubDomains=np.empty((0,4)) if not len(args): Index: /issm/trunk-jpl/src/m/classes/bamgmesh.py =================================================================== --- /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np class bamgmesh(object): @@ -10,20 +10,20 @@ def __init__(self,*args): # {{{ - self.Vertices=numpy.empty((0,3)) - self.Edges=numpy.empty((0,3)) - self.Triangles=numpy.empty((0,0)) - self.Quadrilaterals=numpy.empty((0,0)) - self.IssmEdges=numpy.empty((0,0)) - self.IssmSegments=numpy.empty((0,0)) - self.VerticesOnGeomVertex=numpy.empty((0,0)) - self.VerticesOnGeomEdge=numpy.empty((0,0)) - self.EdgesOnGeomEdge=numpy.empty((0,0)) - self.SubDomains=numpy.empty((0,4)) - self.SubDomainsFromGeom=numpy.empty((0,0)) - self.ElementConnectivity=numpy.empty((0,0)) - self.NodalConnectivity=numpy.empty((0,0)) - self.NodalElementConnectivity=numpy.empty((0,0)) - self.CrackedVertices=numpy.empty((0,0)) - self.CrackedEdges=numpy.empty((0,0)) + self.Vertices=np.empty((0,3)) + self.Edges=np.empty((0,3)) + self.Triangles=np.empty((0,0)) + self.Quadrilaterals=np.empty((0,0)) + self.IssmEdges=np.empty((0,0)) + self.IssmSegments=np.empty((0,0)) + self.VerticesOnGeomVertex=np.empty((0,0)) + self.VerticesOnGeomEdge=np.empty((0,0)) + self.EdgesOnGeomEdge=np.empty((0,0)) + self.SubDomains=np.empty((0,4)) + self.SubDomainsFromGeom=np.empty((0,0)) + self.ElementConnectivity=np.empty((0,0)) + self.NodalConnectivity=np.empty((0,0)) + self.NodalElementConnectivity=np.empty((0,0)) + self.CrackedVertices=np.empty((0,0)) + self.CrackedEdges=np.empty((0,0)) if not len(args): Index: /issm/trunk-jpl/src/m/classes/basalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcings.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/basalforcings.py (revision 21303) @@ -3,5 +3,5 @@ from checkfield import checkfield from WriteData import WriteData -import numpy +import numpy as np class basalforcings(object): @@ -38,10 +38,10 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) print " no basalforcings.groundedice_melting_rate specified: values set as zero" - if numpy.all(numpy.isnan(self.floatingice_melting_rate)): - self.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.floatingice_melting_rate)): + self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices)) print " no basalforcings.floatingice_melting_rate specified: values set as zero" Index: /issm/trunk-jpl/src/m/classes/clusters/vilje.py =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/vilje.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/clusters/vilje.py (revision 21303) @@ -107,9 +107,9 @@ fid.write('#PBS -e %s/%s/%s.errlog \n\n' % (self.executionpath,dirname,modelname)) fid.write('export ISSM_DIR="%s/../"\n' % self.codepath) - fid.write('module load intelcomp/13.0.1\n') - fid.write('module load mpt/2.06\n') - fid.write('module load petsc/3.4.1d\n') - fid.write('module load parmetis/4.0.2\n') - fid.write('module load mumps/4.10.0\n') + fid.write('module load intelcomp/17.0.0\n') + fid.write('module load mpt/2.14\n') + fid.write('module load petsc/3.7.4d\n') + fid.write('module load parmetis/4.0.3\n') + fid.write('module load mumps/5.0.2\n') fid.write('cd %s/%s/\n\n' % (self.executionpath,dirname)) fid.write('mpiexec_mpt -np %i %s/%s %s %s/%s %s\n' % (self.np,self.codepath,executable,str(solution),self.executionpath,dirname,modelname)) Index: /issm/trunk-jpl/src/m/classes/dependent.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dependent.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/dependent.py (revision 21303) @@ -1,4 +1,4 @@ import os.path -import numpy +import numpy as np from pairoptions import pairoptions from fielddisplay import fielddisplay @@ -51,5 +51,5 @@ s+="%s\n" % fielddisplay(self,'type',"type of variable ('vertex' or 'scalar')") - if not numpy.isnan(self.fos_reverse_index): + if not np.isnan(self.fos_reverse_index): s+="%s\n" % fielddisplay(self,'fos_reverse_index',"index for fos_reverse driver of ADOLC") if self.exp: @@ -70,5 +70,5 @@ raise RuntimeError("dependent checkconsistency error: index for segments should be >=0") - if not numpy.isnan(self.fos_reverse_index): + if not np.isnan(self.fos_reverse_index): if not strcmpi(driver,'fos_reverse'): raise TypeError("cannot declare a dependent with a fos_reverse_index when the driver is not fos_reverse!") Index: /issm/trunk-jpl/src/m/classes/esa.py =================================================================== --- /issm/trunk-jpl/src/m/classes/esa.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/esa.py (revision 21303) @@ -2,5 +2,5 @@ from MatlabFuncs import * from model import * -from numpy import * +from np.import * from checkfield import checkfield from WriteData import WriteData Index: /issm/trunk-jpl/src/m/classes/flowequation.py =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import copy from project3d import project3d @@ -105,6 +105,6 @@ md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2]) elif m.strcmp(md.mesh.domaintype(),'3D'): - md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',numpy.arange(0,8+1)) - md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',numpy.arange(0,8+1)) + md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',np.arange(0,8+1)) + md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',np.arange(0,8+1)) else: raise RuntimeError('mesh type not supported yet') @@ -114,5 +114,5 @@ if 'StressbalanceSIAAnalysis' in analyses: if any(self.element_equation==1): - if numpy.any(numpy.logical_and(self.vertex_equation,md.mask.groundedice_levelset)): + if np.any(np.logical_and(self.vertex_equation,md.mask.groundedice_levelset)): print "\n !!! Warning: SIA's model is not consistent on ice shelves !!!\n" Index: /issm/trunk-jpl/src/m/classes/gia.py =================================================================== --- /issm/trunk-jpl/src/m/classes/gia.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/gia.py (revision 21303) @@ -47,6 +47,6 @@ return md - md = checkfield(md,'fieldname','gia.mantle_viscosity','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1],'>',0) - md = checkfield(md,'fieldname','gia.lithosphere_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1],'>',0) + md = checkfield(md,'fieldname','gia.mantle_viscosity','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) + md = checkfield(md,'fieldname','gia.lithosphere_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) md = checkfield(md,'fieldname','gia.cross_section_shape','numel',[1],'values',[1,2]) Index: /issm/trunk-jpl/src/m/classes/groundingline.py =================================================================== --- /issm/trunk-jpl/src/m/classes/groundingline.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/groundingline.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -38,8 +38,8 @@ if not m.strcmp(self.migration,'None'): - if numpy.any(numpy.isnan(md.geometry.bed)): + if np.any(np.isnan(md.geometry.bed)): md.checkmessage("requesting grounding line migration, but bathymetry is absent!") - pos=numpy.nonzero(md.mask.groundedice_levelset>0.)[0] - if any(numpy.abs(md.geometry.base[pos]-md.geometry.bed[pos])>10**-10): + pos=np.nonzero(md.mask.groundedice_levelset>0.)[0] + if any(np.abs(md.geometry.base[pos]-md.geometry.bed[pos])>10**-10): md.checkmessage("base not equal to bed on grounded ice!") if any(md.geometry.bed - md.geometry.base > 10**-9): Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from project3d import project3d from fielddisplay import fielddisplay @@ -134,6 +134,6 @@ # }}} def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.basal_moulin_input)): - self.basal_moulin_input=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.basal_moulin_input)): + self.basal_moulin_input=np.zeros((md.mesh.numberofvertices)) print" no hydrology.basal_moulin_input specified: values set as zero" @@ -166,8 +166,8 @@ md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0.,'numel',[1]) md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0.,'numel',[1]) - md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices]) if self.isefficientlayer==1: md = checkfield(md,'fieldname','hydrology.spcepl_head','Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices,1],'values',[0,1]) + md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices],'values',[0,1]) md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0.,'numel',[1]) md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0.,'numel',[1]) Index: /issm/trunk-jpl/src/m/classes/independent.py =================================================================== --- /issm/trunk-jpl/src/m/classes/independent.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/independent.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from pairoptions import pairoptions from fielddisplay import fielddisplay @@ -17,5 +17,5 @@ self.type = '' self.fos_forward_index = float('NaN') - self.fov_forward_indices = numpy.array([]) + self.fov_forward_indices = np.array([]) self.nods = 0 @@ -34,7 +34,7 @@ s+="%s\n" % fielddisplay(self,'name',"variable name (must match corresponding String)") s+="%s\n" % fielddisplay(self,'type',"type of variable ('vertex' or 'scalar')") - if not numpy.isnan(self.fos_forward_index): + if not np.isnan(self.fos_forward_index): s+="%s\n" % fielddisplay(self,'fos_forward_index',"index for fos_foward driver of ADOLC") - if numpy.any(numpy.logical_not(numpy.isnan(self.fov_forward_indices))): + if np.any(np.logical_not(np.isnan(self.fov_forward_indices))): s+="%s\n" % fielddisplay(self,'fov_forward_indices',"indices for fov_foward driver of ADOLC") @@ -46,5 +46,5 @@ # }}} def checkconsistency(self,md,i,solution,analyses,driver): # {{{ - if not numpy.isnan(self.fos_forward_index): + if not np.isnan(self.fos_forward_index): if not strcmpi(driver,'fos_forward'): raise TypeError("cannot declare an independent with a fos_forward_index when the driver is not fos_forward!") Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from project3d import project3d from fielddisplay import fielddisplay @@ -61,6 +61,13 @@ #Lithostatic pressure by default - self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) - #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1)) + # self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) + #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,)) + + if np.ndim(md.geometry.surface)==2: + print('Reshaping md.geometry.surface for you convenience but you should fix it in you files') + self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface.reshape(-1,)-md.mesh.z) + else: + self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z) + return self #}}} @@ -70,5 +77,5 @@ def checkconsistency(self,md,solution,analyses): # {{{ if 'StressbalanceAnalysis' in analyses: - if not numpy.any(numpy.logical_or(numpy.isnan(md.initialization.vx),numpy.isnan(md.initialization.vy))): + if not np.any(np.logical_or(np.isnan(md.initialization.vx),np.isnan(md.initialization.vy))): md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) @@ -80,6 +87,6 @@ md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) #Triangle with zero velocity - if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ - numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): + if np.any(np.logical_and(np.sum(np.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ + np.sum(np.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): md.checkmessage("at least one triangle has all its vertices with a zero velocity") if 'ThermalAnalysis' in analyses: @@ -98,10 +105,10 @@ if 'HydrologyDCInefficientAnalysis' in analyses: if hasattr(md.hydrology,'hydrologydc'): - md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) if 'HydrologyDCEfficientAnalysis' in analyses: if hasattr(md.hydrology,'hydrologydc'): if md.hydrology.isefficientlayer==1: - md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) - md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) return md @@ -124,7 +131,7 @@ if md.thermal.isenthalpy: tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; - pos = numpy.nonzero(md.initialization.temperature > tpmp)[0] + pos = np.nonzero(md.initialization.temperature > tpmp)[0] enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); - enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,1) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,1) + enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,) WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy'); Index: /issm/trunk-jpl/src/m/classes/inversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/inversion.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from project3d import project3d from fielddisplay import fielddisplay @@ -76,9 +76,9 @@ self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node') self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node') - if not numpy.any(numpy.isnan(self.cost_functions_coefficients)): + if not np.any(np.isnan(self.cost_functions_coefficients)): self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node') - if not numpy.any(numpy.isnan(self.min_parameters)): + if not np.any(np.isnan(self.min_parameters)): self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node') - if not numpy.any(numpy.isnan(self.max_parameters)): + if not np.any(np.isnan(self.max_parameters)): self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node') return self @@ -98,5 +98,5 @@ #maximum number of iteration in the optimization algorithm for #each step - self.maxiter_per_step=20*numpy.ones(self.nsteps) + self.maxiter_per_step=20*np.ones(self.nsteps) #the inversed parameter is updated as follows: @@ -105,5 +105,5 @@ #inversed parameter (10^8 for B, 50 for drag) and can be decreased #after the first iterations - self.gradient_scaling=50*numpy.ones((self.nsteps,1)) + self.gradient_scaling=50*np.ones((self.nsteps,1)) #several responses can be used: @@ -113,5 +113,5 @@ #misfit(1)/misfit(0) < self.step_threshold, we go directly to #the next step - self.step_threshold=.7*numpy.ones(self.nsteps) #30 per cent decrement + self.step_threshold=.7*np.ones(self.nsteps) #30 per cent decrement #cost_function_threshold is a criteria to stop the control methods. @@ -128,6 +128,6 @@ return md - num_controls=numpy.size(md.inversion.control_parameters) - num_costfunc=numpy.size(md.inversion.cost_functions) + num_controls=np.size(md.inversion.control_parameters) + num_costfunc=np.size(md.inversion.cost_functions) md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1]) @@ -185,5 +185,5 @@ #process cost functions - num_cost_functions=numpy.size(self.cost_functions) + num_cost_functions=np.size(self.cost_functions) data=marshallcostfunctions(self.cost_functions) WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') Index: /issm/trunk-jpl/src/m/classes/linearbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/linearbasalforcings.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/linearbasalforcings.py (revision 21303) @@ -2,5 +2,5 @@ from checkfield import checkfield from WriteData import WriteData -import numpy +import numpy as np class linearbasalforcings(object): @@ -51,6 +51,6 @@ def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) print " no basalforcings.groundedice_melting_rate specified: values set as zero" @@ -92,8 +92,8 @@ yts=md.constants.yts - floatingice_melting_rate = numpy.zeros((md.mesh.numberofvertices,1)) - pos=numpy.nonzero(md.geometry.base<=md.basalforcings.deepwater_elevation) + floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) + pos=np.nonzero(md.geometry.base<=md.basalforcings.deepwater_elevation) floatingice_melting_rate[pos]=md.basalforcings.deepwater_melting_rate - pos=numpy.nonzero(numpy.logical_and(md.geometry.base>md.basalforcings.deepwater_elevation,md.geometry.basemd.basalforcings.deepwater_elevation,md.geometry.base 0, grounding line position if = 0, floating ice if < 0')) string="%s\n%s"%(string,fielddisplay(self,'ice_levelset','presence of ice if < 0, icefront position if = 0, no ice if > 0')) string="%s\n%s"%(string,fielddisplay(self,'ocean_levelset','is the vertex on the ocean ? yes if = 1, no if = 0')) string="%s\n%s"%(string,fielddisplay(self,'land_levelset','is the vertex on the land ? yes if = 1, no if = 0')) - return string + # }}} def loadobj(self): # {{{ # This def is directly called by matlab when a model object is # loaded. Update old properties here - #2014 February 5th if numel(self.ice_levelset)>1 and all(self.ice_levelset>=0): @@ -41,13 +40,15 @@ return self # }}} + def setdefaultparameters(self): # {{{ return self + # }}} + def checkconsistency(self,md,solution,analyses): # {{{ - - md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices, 1]) - md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices, 1]) - md = checkfield(md,'fieldname','mask.ocean_levelset','size',[md.mesh.numberofvertices, 1]) - md = checkfield(md,'fieldname','mask.land_levelset','size',[md.mesh.numberofvertices, 1]) + md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mask.ocean_levelset','size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','mask.land_levelset','size',[md.mesh.numberofvertices]) isice=(md.mask.ice_levelset<=0) if sum(isice)==0: @@ -56,12 +57,14 @@ if max(md.mask.ice_levelset)<0: print('no ice front provided') - - elements=md.mesh.elements-1; elements=elements.astype(numpy.int32, copy=False); - icefront=numpy.sum(md.mask.ice_levelset[elements]==0,axis=1) + + elements=md.mesh.elements-1; elements=elements.astype(np.int32, copy=False); + icefront=np.sum(md.mask.ice_levelset[elements]==0,axis=1) if (max(icefront)==3 & m.strcmp(md.mesh.elementtype(),'Tria')) or (max(icefront==6) & m.strcmp(md.mesh.elementtype(),'Penta')): raise RuntimeError('At least one element has all nodes on ice front, change md.mask.ice_levelset to fix it') - + return md + # }}} + def extrude(self,md): # {{{ self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node') @@ -71,4 +74,5 @@ return self # }}} + def mask(*args): # {{{ if not len(args): @@ -77,5 +81,7 @@ raise RuntimeError('constructor not supported') return self + # }}} + def marshall(self,prefix,md,fid): # {{{ WriteData(fid,prefix,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1) @@ -85,9 +91,8 @@ # }}} def savemodeljs(self,fid,modelname): # {{{ - writejs1Darray(fid,[modelname, '.mask.groundedice_levelset'],self.groundedice_levelset) writejs1Darray(fid,[modelname, '.mask.ice_levelset'],self.ice_levelset) writejs1Darray(fid,[modelname, '.mask.ocean_levelset'],self.ocean_levelset) writejs1Darray(fid,[modelname, '.mask.land_levelset'],self.land_levelset) + # }}} - # }}} Index: /issm/trunk-jpl/src/m/classes/mesh2d.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh2d.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/mesh2d.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -85,7 +85,7 @@ md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1)) + md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1)) md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3]) - if numpy.any(numpy.logical_not(m.ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): + if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): md.checkmessage("orphan nodes have been found. Check the mesh outline") md = checkfield(md,'fieldname','mesh.numberofelements','>',0) @@ -112,5 +112,5 @@ WriteData(fid,prefix,'object',self,'class','mesh','fieldname','x','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'class','mesh','fieldname','y','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'name','md.mesh.z','data',numpy.zeros(self.numberofvertices),'format','DoubleMat','mattype',1); + WriteData(fid,prefix,'name','md.mesh.z','data',np.zeros(self.numberofvertices),'format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2) WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofelements','format','Integer') Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dprisms.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import * @@ -108,7 +108,7 @@ md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1)) + md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1)) md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,6]) - if numpy.any(numpy.logical_not(m.ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): + if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): md.checkmessage("orphan nodes have been found. Check the mesh3dprisms outline") md = checkfield(md,'fieldname','mesh.numberoflayers','>=',0) Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 21303) @@ -1,5 +1,5 @@ from MatlabFuncs import * from model import * -from numpy import * +from np.import * from fielddisplay import fielddisplay from checkfield import checkfield Index: /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/mismipbasalforcings.py (revision 21303) @@ -3,5 +3,5 @@ from checkfield import checkfield from WriteData import WriteData -import numpy +import numpy as np class mismipbasalforcings(object): @@ -21,6 +21,6 @@ self.geothermalflux = float('NaN') - if numpy.all(numpy.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate=numpy.zeros(md.mesh.numberofvertices) + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate=np.zeros(md.mesh.numberofvertices) print ' no basalforcings.groundedice_melting_rate specified: values set as zero' @@ -81,6 +81,6 @@ print 'WARNING: value of yts for MISMIP+ runs different from ISSM default!' - floatingice_melting_rate = numpy.zeros((md.mesh.numberofvertices,1)) - floatingice_melting_rate = md.basalforcings.meltrate_factor*numpy.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*numpy.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0) + floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) + floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*np.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0) WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 21303) @@ -1,4 +1,4 @@ #module imports {{{ -import numpy +import numpy as np import copy import sys @@ -231,32 +231,32 @@ #get elements that are inside area flag_elem=FlagElements(md1,area) - if not numpy.any(flag_elem): + if not np.any(flag_elem): raise RuntimeError("extracted model is empty") #kick out all elements with 3 dirichlets - spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] - spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1 - flag=numpy.ones(md1.mesh.numberofvertices) + spc_elem=np.nonzero(np.logical_not(flag_elem))[0] + spc_node=np.unique(md1.mesh.elements[spc_elem,:])-1 + flag=np.ones(md1.mesh.numberofvertices) flag[spc_node]=0 - pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0] + pos=np.nonzero(np.logical_not(np.sum(flag[md1.mesh.elements-1],axis=1)))[0] flag_elem[pos]=0 #extracted elements and nodes lists - pos_elem=numpy.nonzero(flag_elem)[0] - pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1 + pos_elem=np.nonzero(flag_elem)[0] + pos_node=np.unique(md1.mesh.elements[pos_elem,:])-1 #keep track of some fields numberofvertices1=md1.mesh.numberofvertices numberofelements1=md1.mesh.numberofelements - numberofvertices2=numpy.size(pos_node) - numberofelements2=numpy.size(pos_elem) - flag_node=numpy.zeros(numberofvertices1) + numberofvertices2=np.size(pos_node) + numberofelements2=np.size(pos_elem) + flag_node=np.zeros(numberofvertices1) flag_node[pos_node]=1 #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) - Pelem=numpy.zeros(numberofelements1,int) - Pelem[pos_elem]=numpy.arange(1,numberofelements2+1) - Pnode=numpy.zeros(numberofvertices1,int) - Pnode[pos_node]=numpy.arange(1,numberofvertices2+1) + Pelem=np.zeros(numberofelements1,int) + Pelem[pos_elem]=np.arange(1,numberofelements2+1) + Pnode=np.zeros(numberofvertices1,int) + Pnode[pos_node]=np.arange(1,numberofvertices2+1) #renumber the elements (some node won't exist anymore) @@ -283,5 +283,5 @@ #get field field=getattr(md1,fieldi) - fieldsize=numpy.shape(field) + fieldsize=np.shape(field) if hasattr(field,'__dict__') and not m.ismember(fieldi,['results'])[0]: #recursive call object_fields=vars(field) @@ -289,5 +289,5 @@ #get field field=getattr(getattr(md1,fieldi),fieldj) - fieldsize=numpy.shape(field) + fieldsize=np.shape(field) if len(fieldsize): #size = number of nodes * n @@ -295,5 +295,5 @@ setattr(getattr(md2,fieldi),fieldj,field[pos_node]) elif fieldsize[0]==numberofvertices1+1: - setattr(getattr(md2,fieldi),fieldj,numpy.vstack((field[pos_node],field[-1,:]))) + setattr(getattr(md2,fieldi),fieldj,np.vstack((field[pos_node],field[-1,:]))) #size = number of elements * n elif fieldsize[0]==numberofelements1: @@ -305,5 +305,5 @@ setattr(md2,fieldi,field[pos_node]) elif fieldsize[0]==numberofvertices1+1: - setattr(md2,fieldi,numpy.hstack((field[pos_node],field[-1,:]))) + setattr(md2,fieldi,np.hstack((field[pos_node],field[-1,:]))) #size = number of elements * n elif fieldsize[0]==numberofelements1: @@ -320,28 +320,28 @@ if md1.mesh.__class__.__name__=='mesh3dprisms': md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node] - pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0] - md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1] + pos=np.where(~np.isnan(md2.mesh.uppervertex))[0] + md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1] md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node] - pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0] - md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1] + pos=np.where(~np.isnan(md2.mesh.lowervertex))[0] + md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1] md2.mesh.upperelements=md1.mesh.upperelements[pos_elem] - pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0] - md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1] + pos=np.where(~np.isnan(md2.mesh.upperelements))[0] + md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1] md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem] - pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0] - md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1] + pos=np.where(~np.isnan(md2.mesh.lowerelements))[0] + md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1] #Initial 2d mesh if md1.mesh.__class__.__name__=='mesh3dprisms': - flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)] - pos_elem_2d=numpy.nonzero(flag_elem_2d)[0] - flag_node_2d=flag_node[numpy.arange(0,md1.mesh.numberofvertices2d)] - pos_node_2d=numpy.nonzero(flag_node_2d)[0] - - md2.mesh.numberofelements2d=numpy.size(pos_elem_2d) - md2.mesh.numberofvertices2d=numpy.size(pos_node_2d) + flag_elem_2d=flag_elem[np.arange(0,md1.mesh.numberofelements2d)] + pos_elem_2d=np.nonzero(flag_elem_2d)[0] + flag_node_2d=flag_node[np.arange(0,md1.mesh.numberofvertices2d)] + pos_node_2d=np.nonzero(flag_node_2d)[0] + + md2.mesh.numberofelements2d=np.size(pos_elem_2d) + md2.mesh.numberofvertices2d=np.size(pos_node_2d) md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:] md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1] @@ -354,7 +354,7 @@ #Edges if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'): - if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some numpy.nans... + if np.ndim(md2.mesh.edges)>1 and np.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some np.nans... #renumber first two columns - pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0] + pos=np.nonzero(md2.mesh.edges[:,3]!=-1)[0] md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1] md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1]-1] @@ -362,12 +362,12 @@ md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1] #remove edges when the 2 vertices are not in the domain. - md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:] + md2.mesh.edges=md2.mesh.edges[np.nonzero(np.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:] #Replace all zeros by -1 in the last two columns - pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0] + pos=np.nonzero(md2.mesh.edges[:,2]==0)[0] md2.mesh.edges[pos,2]=-1 - pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0] + pos=np.nonzero(md2.mesh.edges[:,3]==0)[0] md2.mesh.edges[pos,3]=-1 #Invert -1 on the third column with last column (Also invert first two columns!!) - pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0] + pos=np.nonzero(md2.mesh.edges[:,2]==-1)[0] md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3] md2.mesh.edges[pos,3]=-1 @@ -376,16 +376,16 @@ md2.mesh.edges[pos,0]=values #Finally remove edges that do not belong to any element - pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0] - md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0) + pos=np.nonzero(np.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0] + md2.mesh.edges=np.delete(md2.mesh.edges,pos,axis=0) #Penalties - if numpy.any(numpy.logical_not(numpy.isnan(md2.stressbalance.vertex_pairing))): - for i in xrange(numpy.size(md1.stressbalance.vertex_pairing,axis=0)): + if np.any(np.logical_not(np.isnan(md2.stressbalance.vertex_pairing))): + for i in xrange(np.size(md1.stressbalance.vertex_pairing,axis=0)): md2.stressbalance.vertex_pairing[i,:]=Pnode[md1.stressbalance.vertex_pairing[i,:]] - md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[numpy.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:] - if numpy.any(numpy.logical_not(numpy.isnan(md2.masstransport.vertex_pairing))): - for i in xrange(numpy.size(md1.masstransport.vertex_pairing,axis=0)): + md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[np.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:] + if np.any(np.logical_not(np.isnan(md2.masstransport.vertex_pairing))): + for i in xrange(np.size(md1.masstransport.vertex_pairing,axis=0)): md2.masstransport.vertex_pairing[i,:]=Pnode[md1.masstransport.vertex_pairing[i,:]] - md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[numpy.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:] + md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[np.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:] #recreate segments @@ -394,5 +394,5 @@ md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)[0] md2.mesh.segments=contourenvelope(md2) - md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool) + md2.mesh.vertexonboundary=np.zeros(numberofvertices2,bool) md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True else: @@ -401,7 +401,7 @@ md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)[0] segments=contourenvelope(md2) - md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) + md2.mesh.vertexonboundary=np.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) md2.mesh.vertexonboundary[segments[:,0:2]-1]=True - md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) + md2.mesh.vertexonboundary=np.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) #Then do it for 3d as usual md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)[0] @@ -410,21 +410,21 @@ #Boundary conditions: Dirichlets on new boundary #Catch the elements that have not been extracted - orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] - orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1 + orphans_elem=np.nonzero(np.logical_not(flag_elem))[0] + orphans_node=np.unique(md1.mesh.elements[orphans_elem,:])-1 #Figure out which node are on the boundary between md2 and md1 - nodestoflag1=numpy.intersect1d(orphans_node,pos_node) + nodestoflag1=np.intersect1d(orphans_node,pos_node) nodestoflag2=Pnode[nodestoflag1].astype(int)-1 - if numpy.size(md1.stressbalance.spcvx)>1 and numpy.size(md1.stressbalance.spcvy)>2 and numpy.size(md1.stressbalance.spcvz)>2: - if numpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1: + if np.size(md1.stressbalance.spcvx)>1 and np.size(md1.stressbalance.spcvy)>2 and np.size(md1.stressbalance.spcvz)>2: + if np.size(md1.inversion.vx_obs)>1 and np.size(md1.inversion.vy_obs)>1: md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2] else: - md2.stressbalance.spcvx[nodestoflag2]=numpy.nan - md2.stressbalance.spcvy[nodestoflag2]=numpy.nan + md2.stressbalance.spcvx[nodestoflag2]=np.nan + md2.stressbalance.spcvy[nodestoflag2]=np.nan print "\n!! extract warning: spc values should be checked !!\n\n" #put 0 for vz md2.stressbalance.spcvz[nodestoflag2]=0 - if numpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))): - md2.thermal.spctemperature[nodestoflag2,0]=1 + if np.any(np.logical_not(np.isnan(md1.thermal.spctemperature))): + md2.thermal.spctemperature[nodestoflag2]=1 #Results fields @@ -441,7 +441,7 @@ #get subfields for solutionsubfield,subfield in fieldi.__dict__.iteritems(): - if numpy.size(subfield)==numberofvertices1: + if np.size(subfield)==numberofvertices1: setattr(fieldr,solutionsubfield,subfield[pos_node]) - elif numpy.size(subfield)==numberofelements1: + elif np.size(subfield)==numberofelements1: setattr(fieldr,solutionsubfield,subfield[pos_elem]) else: @@ -455,7 +455,7 @@ #get subfields for solutionsubfield,subfield in field.__dict__.iteritems(): - if numpy.size(subfield)==numberofvertices1: + if np.size(subfield)==numberofvertices1: setattr(fieldr,solutionsubfield,subfield[pos_node]) - elif numpy.size(subfield)==numberofelements1: + elif np.size(subfield)==numberofelements1: setattr(fieldr,solutionsubfield,subfield[pos_elem]) else: @@ -510,5 +510,5 @@ raise TypeError("extrusionexponent must be >=0") numlayers=args[0] - extrusionlist=(numpy.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1] + extrusionlist=(np.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1] elif len(args)==3: #two polynomial laws @@ -520,7 +520,7 @@ raise TypeError("lower and upper extrusionexponents must be >=0") - lowerextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2. - upperextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2. - extrusionlist=numpy.unique(numpy.concatenate((lowerextrusionlist,1.-upperextrusionlist))) + lowerextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2. + upperextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2. + extrusionlist=np.unique(np.concatenate((lowerextrusionlist,1.-upperextrusionlist))) if numlayers<2: @@ -550,7 +550,7 @@ md.mesh.extractedelements = mesh2d.extractedelements - x3d=numpy.empty((0)) - y3d=numpy.empty((0)) - z3d=numpy.empty((0)) #the lower node is on the bed + x3d=np.empty((0)) + y3d=np.empty((0)) + z3d=np.empty((0)) #the lower node is on the bed thickness3d=md.geometry.thickness #thickness and bed for these nodes bed3d=md.geometry.base @@ -558,29 +558,29 @@ #Create the new layers for i in xrange(numlayers): - x3d=numpy.concatenate((x3d,md.mesh.x)) - y3d=numpy.concatenate((y3d,md.mesh.y)) + x3d=np.concatenate((x3d,md.mesh.x)) + y3d=np.concatenate((y3d,md.mesh.y)) #nodes are distributed between bed and surface accordingly to the given exponent - z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1))) - number_nodes3d=numpy.size(x3d) #number of 3d nodes for the non extruded part of the mesh + z3d=np.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1))) + number_nodes3d=np.size(x3d) #number of 3d nodes for the non extruded part of the mesh #Extrude elements - elements3d=numpy.empty((0,6),int) + elements3d=np.empty((0,6),int) for i in xrange(numlayers-1): - elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices)))) #Create the elements of the 3d mesh for the non extruded part - number_el3d=numpy.size(elements3d,axis=0) #number of 3d nodes for the non extruded part of the mesh + elements3d=np.vstack((elements3d,np.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices)))) #Create the elements of the 3d mesh for the non extruded part + number_el3d=np.size(elements3d,axis=0) #number of 3d nodes for the non extruded part of the mesh #Keep a trace of lower and upper nodes - lowervertex=-1*numpy.ones(number_nodes3d,int) - uppervertex=-1*numpy.ones(number_nodes3d,int) - lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1) - uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1) + lowervertex=np.nan*np.ones(number_nodes3d,int) + uppervertex=np.nan*np.ones(number_nodes3d,int) + lowervertex[md.mesh.numberofvertices:]=np.arange(1,(numlayers-1)*md.mesh.numberofvertices+1) + uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=np.arange(md.mesh.numberofvertices+1,number_nodes3d+1) md.mesh.lowervertex=lowervertex md.mesh.uppervertex=uppervertex #same for lower and upper elements - lowerelements=-1*numpy.ones(number_el3d,int) - upperelements=-1*numpy.ones(number_el3d,int) - lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1) - upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1) + lowerelements=np.nan*np.ones(number_el3d,int) + upperelements=np.nan*np.ones(number_el3d,int) + lowerelements[md.mesh.numberofelements:]=np.arange(1,(numlayers-2)*md.mesh.numberofelements+1) + upperelements[:(numlayers-2)*md.mesh.numberofelements]=np.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1) md.mesh.lowerelements=lowerelements md.mesh.upperelements=upperelements @@ -605,6 +605,6 @@ #bedinfo and surface info - md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1) - md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers) + md.mesh.vertexonbase=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1) + md.mesh.vertexonsurface=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers) md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node') @@ -630,11 +630,11 @@ #connectivity - md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1)) - md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1 - if not numpy.isnan(md.mesh.elementconnectivity).all(): + md.mesh.elementconnectivity=np.tile(md.mesh.elementconnectivity,(numlayers-1,1)) + md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1 + if not np.isnan(md.mesh.elementconnectivity).all(): for i in xrange(1,numlayers-1): md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \ =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d - md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0 + md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity<0)]=0 md.materials.extrude(md) @@ -674,32 +674,32 @@ #observations - if not numpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) - if not numpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) - if not numpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) - if not numpy.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) - if isinstance(md.inversion.min_parameters,numpy.ndarray): + if not np.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) + if isinstance(md.inversion.min_parameters,np.ndarray): if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) - if isinstance(md.inversion.max_parameters,numpy.ndarray): + if isinstance(md.inversion.max_parameters,np.ndarray): if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers) - if not numpy.isnan(md.smb.mass_balance).all(): + if not np.isnan(md.smb.mass_balance).all(): md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) #results - if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) - if not numpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) - if not numpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) - if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) - if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) - if not numpy.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1) - if not numpy.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) - if not numpy.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1) - if not numpy.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1) + if not np.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) + if not np.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) + if not np.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) + if not np.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) + if not np.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) + if not np.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1) + if not np.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) + if not np.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1) + if not np.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1) #gia - if not numpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) - if not numpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) + if not np.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) + if not np.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) #elementstype - if not numpy.isnan(md.flowequation.element_equation).all(): + if not np.isnan(md.flowequation.element_equation).all(): md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1) md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1) @@ -725,5 +725,5 @@ md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers) md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers) - if not numpy.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) + if not np.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1) @@ -752,5 +752,5 @@ md.geometry.thickness=project2d(md,md.geometry.thickness,1) md.geometry.base=project2d(md,md.geometry.base,1) - if isinstance(md.geometry.bed,numpy.ndarray): + if isinstance(md.geometry.bed,np.ndarray): md.geometry.bed=project2d(md,md.geometry.bed,1) md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1) @@ -758,7 +758,7 @@ #lat long - if isinstance(md.mesh.lat,numpy.ndarray): + if isinstance(md.mesh.lat,np.ndarray): if md.mesh.lat.size==md.mesh.numberofvertices: md.mesh.lat=project2d(md,md.mesh.lat,1) - if isinstance(md.mesh.long,numpy.ndarray): + if isinstance(md.mesh.long,np.ndarray): if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1) @@ -770,6 +770,6 @@ mesh.numberofelements=md.mesh.numberofelements2d mesh.elements=md.mesh.elements2d - if not numpy.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) - if not numpy.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) + if not np.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) + if not np.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) md.mesh=mesh Index: /issm/trunk-jpl/src/m/classes/outputdefinition.py =================================================================== --- /issm/trunk-jpl/src/m/classes/outputdefinition.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/outputdefinition.py (revision 21303) @@ -2,5 +2,5 @@ from checkfield import checkfield from WriteData import WriteData -import numpy as np +import numpy as np class outputdefinition(object): Index: /issm/trunk-jpl/src/m/classes/qmu.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/qmu.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from project3d import project3d from collections import OrderedDict @@ -105,5 +105,5 @@ # }}} def extrude(self,md): # {{{ - self.partition=project3d(md,'vector',numpy.transpose(self.partition),'type','node') + self.partition=project3d(md,'vector',np.transpose(self.partition),'type','node') return self #}}} @@ -120,5 +120,5 @@ md.checkmessage("concurrency should be set to 1 when running dakota in library mode") if md.qmu.partition: - if not numpy.size(md.qmu.partition)==md.mesh.numberofvertices: + if not np.size(md.qmu.partition)==md.mesh.numberofvertices: md.checkmessage("user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1") if not min(md.qmu.partition)==0: Index: /issm/trunk-jpl/src/m/classes/results.py =================================================================== --- /issm/trunk-jpl/src/m/classes/results.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/results.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from pairoptions import pairoptions from fielddisplay import fielddisplay Index: /issm/trunk-jpl/src/m/classes/rifts.py =================================================================== --- /issm/trunk-jpl/src/m/classes/rifts.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/rifts.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -33,5 +33,5 @@ #}}} def checkconsistency(self,md,solution,analyses): # {{{ - if (not self.riftstruct) or numpy.any(isnans(self.riftstruct)): + if (not self.riftstruct) or np.any(isnans(self.riftstruct)): numrifts=0 else: @@ -43,5 +43,5 @@ if not isinstance(self.riftstruct,list): md.checkmessage("rifts.riftstruct should be a structure!") - if numpy.any(md.mesh.segmentmarkers>=2): + if np.any(md.mesh.segmentmarkers>=2): #We have segments with rift markers, but no rift structure! md.checkmessage("model should be processed for rifts (run meshprocessrifts)!") @@ -49,5 +49,5 @@ md = checkfield(md,'fieldname',"rifts.riftstruct[%d]['fill']" % i,'values',['Water','Air','Ice','Melange',0,1,2,3]) else: - if self.riftstruct and numpy.any(numpy.logical_not(isnans(self.riftstruct))): + if self.riftstruct and np.any(np.logical_not(isnans(self.riftstruct))): md.checkmessage("riftstruct should be NaN since numrifts is 0!") @@ -57,5 +57,5 @@ #Process rift info - if (not self.riftstruct) or numpy.any(isnans(self.riftstruct)): + if (not self.riftstruct) or np.any(isnans(self.riftstruct)): numrifts=0 else: @@ -64,5 +64,5 @@ numpairs=0 for rift in self.riftstruct: - numpairs+=numpy.size(rift['penaltypairs'],axis=0) + numpairs+=np.size(rift['penaltypairs'],axis=0) # Convert strings in riftstruct to hard coded numbers @@ -76,8 +76,8 @@ # 2 for nodes + 2 for elements+ 2 for normals + 1 for length + 1 for fill + 1 for friction + 1 for fraction + 1 for fractionincrement + 1 for state. - data=numpy.zeros((numpairs,12)) + data=np.zeros((numpairs,12)) count=0 for rift in self.riftstruct: - numpairsforthisrift=numpy.size(rift['penaltypairs'],0) + numpairsforthisrift=np.size(rift['penaltypairs'],0) data[count:count+numpairsforthisrift,0:7]=rift['penaltypairs'] data[count:count+numpairsforthisrift,7]=rift['fill'] Index: /issm/trunk-jpl/src/m/classes/slr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/slr.py (revision 21303) @@ -2,5 +2,5 @@ from MatlabFuncs import * from model import * -from numpy import * +import numpy as np from checkfield import checkfield from WriteData import WriteData @@ -15,6 +15,6 @@ def __init__(self): # {{{ - self.deltathickness = NaN - self.sealevel = NaN + self.deltathickness = float('NaN') + self.sealevel = float('NaN') self.maxiter = 0 self.reltol = 0 @@ -60,5 +60,5 @@ #Convergence criterion: absolute, relative and residual - self.reltol=NaN #default + self.reltol=float('NaN') #default self.abstol=0.001 #1 mm of sea level rise @@ -95,6 +95,6 @@ return md - md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements,1]) - md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1]) + md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) + md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1) md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1) Index: /issm/trunk-jpl/src/m/classes/steadystate.py =================================================================== --- /issm/trunk-jpl/src/m/classes/steadystate.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/steadystate.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from fielddisplay import fielddisplay from checkfield import checkfield @@ -54,5 +54,5 @@ md.checkmessage("for a steadystate computation, timestepping.time_step must be zero.") - if numpy.isnan(md.stressbalance.reltol): + if np.isnan(md.stressbalance.reltol): md.checkmessage("for a steadystate computation, stressbalance.reltol (relative convergence criterion) must be defined!") Index: /issm/trunk-jpl/src/m/classes/stressbalance.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import sys import copy @@ -149,24 +149,24 @@ #singular solution # if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2), - if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)),numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)))): + if not np.any(np.logical_and(np.logical_not(np.isnan(md.stressbalance.spcvx)),np.logical_not(np.isnan(md.stressbalance.spcvy)))): print "\n !!! Warning: no spc applied, model might not be well posed if no basal friction is applied, check for solution crash\n" #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES # if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6), - if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=6)): + if np.any(np.logical_and(np.sum(np.isnan(md.stressbalance.referential),axis=1)!=0,np.sum(np.isnan(md.stressbalance.referential),axis=1)!=6)): md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values") #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL # if any(sum(isnan(md.stressbalance.referential),2)==0), - if numpy.any(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)==0): - pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)) if item==0] -# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply -# if numpy.any(numpy.abs(numpy.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon): + if np.any(np.sum(np.isnan(md.stressbalance.referential),axis=1)==0): + pos=[i for i,item in enumerate(np.sum(np.isnan(md.stressbalance.referential),axis=1)) if item==0] +# np.inner (and np.dot) calculate all the dot product permutations, resulting in a full matrix multiply +# if np.any(np.abs(np.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon): # md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal") for item in md.stressbalance.referential[pos,:]: - if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon: + if np.abs(np.inner(item[0:2],item[3:5]))>sys.float_info.epsilon: md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal") #CHECK THAT NO rotation specified for FS Grounded ice at base if m.strcmp(md.mesh.domaintype(),'3D') and md.flowequation.isFS: - pos=numpy.nonzero(numpy.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase)) - if numpy.any(numpy.logical_not(numpy.isnan(md.stressbalance.referential[pos,:]))): + pos=np.nonzero(np.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase)) + if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos,:]))): md.checkmessage("no referential should be specified for basal vertices of grounded ice") @@ -195,5 +195,5 @@ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1) - if isinstance(self.loadingforce, (list, tuple, numpy.ndarray)): + if isinstance(self.loadingforce, (list, tuple, np.ndarray)): lx=self.loadingforce[:,0]; ly=self.loadingforce[:,1]; Index: /issm/trunk-jpl/src/m/classes/taoinversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/taoinversion.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/taoinversion.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from project3d import project3d from WriteData import WriteData @@ -125,6 +125,6 @@ - num_controls= numpy.numel(md.inversion.control_parameters) - num_costfunc= numpy.size(md.inversion.cost_functions,2) + num_controls= np.numel(md.inversion.control_parameters) + num_costfunc= np.size(md.inversion.cost_functions,2) md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1]) @@ -155,10 +155,10 @@ if solution=='BalancethicknessSolution': - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) elif solution=='BalancethicknessSoftSolution': - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) else: - md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) def marshall(self, md, fid): @@ -188,10 +188,10 @@ #process control parameters - num_control_parameters = numpy.numel(self.control_parameters) + num_control_parameters = np.numel(self.control_parameters) WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray') WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer') #process cost functions - num_cost_functions = numpy.size(self.cost_functions,2) + num_cost_functions = np.size(self.cost_functions,2) data= marshallcostfunctions(self.cost_functions) WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') Index: /issm/trunk-jpl/src/m/classes/thermal.py =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.py (revision 21302) +++ /issm/trunk-jpl/src/m/classes/thermal.py (revision 21303) @@ -1,8 +1,7 @@ -import numpy +import numpy as np from project3d import project3d from fielddisplay import fielddisplay from checkfield import checkfield from WriteData import WriteData -import MatlabFuncs as m class thermal(object): @@ -44,8 +43,8 @@ #}}} def extrude(self,md): # {{{ - self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan) - if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: - self.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) - pos=numpy.nonzero(md.mesh.vertexonsurface)[0] + self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',np.nan) + if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: + self.spctemperature=float('NaN')*np.ones((md.mesh.numberofvertices)) + pos=np.where(md.mesh.vertexonsurface)[0] self.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface return self @@ -95,15 +94,22 @@ md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2]) md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1) + md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) + if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3: - pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))) - replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) - md = checkfield(md,'fieldname','thermal.spctemperature[numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices,:])))]','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate[pos],'message',"spctemperature should be below the adjusted melting point") + pos=np.where(~np.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])) + try: + spccol=np.size(md.thermal.spctemperature,1) + except IndexError: + spccol=1 + replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)) + control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate + md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1]) md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); if(md.thermal.isenthalpy): - if numpy.isnan(md.stressbalance.reltol): + if np.isnan(md.stressbalance.reltol): md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!") md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero"); - md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) + return md Index: /issm/trunk-jpl/src/m/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 21302) +++ /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os from pairoptions import pairoptions @@ -44,5 +44,5 @@ if isinstance(field,(bool,int,long,float)): - field=numpy.array([field]) + field=np.array([field]) #check empty @@ -55,21 +55,27 @@ if options.exist('size'): fieldsize=options.getfieldvalue('size') - if len(fieldsize) == 1: - if numpy.isnan(fieldsize[0]): + if len(fieldsize) == 1: + if np.isnan(fieldsize[0]): pass - elif not numpy.size(field,0)==fieldsize[0]: - md = md.checkmessage(options.getfieldvalue('message',\ - "field '%s' size should be %d" % (fieldname,fieldsize[0]))) + elif np.ndim(field)==1: + if not np.size(field)==fieldsize[0]: + md = md.checkmessage(options.getfieldvalue('message',"field {} size should be {}".format(fieldname,fieldsize[0]))) + else: + try: + exec("md.{}=field[:,0]".format(fieldname)) + print('{} had a bad dimension, we fixed it but you should check it'.format(fieldname)) + except IndexError: + md = md.checkmessage(options.getfieldvalue('message',"field {} should have {} dimension".format(fieldname,len(fieldsize)))) elif len(fieldsize) == 2: - if numpy.isnan(fieldsize[0]): - if not numpy.size(field,1)==fieldsize[1]: + if np.isnan(fieldsize[0]): + if not np.size(field,1)==fieldsize[1]: md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) - elif numpy.isnan(fieldsize[1]): - if not numpy.size(field,0)==fieldsize[0]: + elif np.isnan(fieldsize[1]): + if not np.size(field,0)==fieldsize[0]: md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) else: - if (not numpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]): + if (not np.size(field,0)==fieldsize[0]) or (not np.size(field,1)==fieldsize[1]): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) @@ -78,5 +84,5 @@ if options.exist('numel'): fieldnumel=options.getfieldvalue('numel') - if numpy.size(field) not in fieldnumel: + if np.size(field) not in fieldnumel: if len(fieldnumel)==1: md = md.checkmessage(options.getfieldvalue('message',\ @@ -91,5 +97,5 @@ #check NaN if options.getfieldvalue('NaN',0): - if numpy.any(numpy.isnan(field)): + if np.any(np.isnan(field)): md = md.checkmessage(options.getfieldvalue('message',\ "NaN values found in field '%s'" % fieldname)) @@ -97,5 +103,5 @@ #check Inf if options.getfieldvalue('Inf',0): - if numpy.any(numpy.isinf(field)): + if np.any(np.isinf(field)): md = md.checkmessage(options.getfieldvalue('message',\ "Inf values found in field '%s'" % fieldname)) @@ -124,10 +130,10 @@ if options.exist('>='): lowerbound=options.getfieldvalue('>=') - if numpy.any(field'): lowerbound=options.getfieldvalue('>') - if numpy.any(field<=lowerbound): + if np.any(field<=lowerbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values above %d" % (fieldname,lowerbound))) @@ -136,10 +142,10 @@ if options.exist('<='): upperbound=options.getfieldvalue('<=') - if numpy.any(field>upperbound): + if np.any(field>upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) if options.exist('<'): upperbound=options.getfieldvalue('<') - if numpy.any(field>=upperbound): + if np.any(field>=upperbound): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have values below %d" % (fieldname,upperbound))) @@ -158,10 +164,10 @@ #Check forcings (size and times) if options.getfieldvalue('timeseries',0): - if numpy.size(field,0)==md.mesh.numberofvertices: - if numpy.ndim(field)>1 and not numpy.size(field,1)==1: + if np.size(field,0)==md.mesh.numberofvertices: + if np.ndim(field)>1 and not np.size(field,1)==1: md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) - elif numpy.size(field,0)==md.mesh.numberofvertices+1 or numpy.size(field,0)==2: - if not all(field[-1,:]==numpy.sort(field[-1,:])): + elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==2: + if not all(field[-1,:]==np.sort(field[-1,:])): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' columns should be sorted chronologically" % fieldname)) @@ -175,6 +181,6 @@ #Check single value forcings (size and times) if options.getfieldvalue('singletimeseries',0): - if numpy.size(field,0)==2: - if not all(field[-1,:]==numpy.sort(field[-1,:])): + if np.size(field,0)==2: + if not all(field[-1,:]==np.sort(field[-1,:])): md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' columns should be sorted chronologically" % fieldname)) Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 21302) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 21303) @@ -1,4 +1,4 @@ from netCDF4 import Dataset, stringtochar -import numpy as np +import numpy as np import time import collections @@ -42,5 +42,5 @@ typelist=[bool,str,unicode,int,float,complex, collections.OrderedDict, - numpy.int64,numpy.ndarray,numpy.float64] + np.int64,np.ndarray,np.float64] groups=dict.keys(md.__dict__) for group in groups: @@ -121,8 +121,8 @@ TypeDict = {float:'f8', 'float64':'f8', - numpy.float64:'f8', + np.float64:'f8', int:'i8', 'int64':'i8', - numpy.int64:'i8', + np.int64:'i8', str:str, dict:str} Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 21302) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os import model @@ -40,18 +40,18 @@ os.mkdir(filename) - IsEnveloppe=numpy.where(model.mesh.vertexonbase | model.mesh.vertexonsurface) + IsEnveloppe=np.where(model.mesh.vertexonbase | model.mesh.vertexonsurface) #get the element related variables if 'z' in dict.keys(model.mesh.__dict__): - points=numpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) - num_of_elt=numpy.size(numpy.isnan(model.mesh.lowerelements))+numpy.size(numpy.isnan(model.mesh.upperelements)) - low_elt_num=numpy.size(numpy.isnan(model.mesh.lowerelements)) - top_elt_num=numpy.size(numpy.isnan(model.mesh.upperelements)) + points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) + num_of_elt=np.size(np.isnan(model.mesh.lowerelements))+np.size(np.isnan(model.mesh.upperelements)) + low_elt_num=np.size(np.isnan(model.mesh.lowerelements)) + top_elt_num=np.size(np.isnan(model.mesh.upperelements)) else: - points=numpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x)))) - num_of_elt=numpy.shape(model.mesh.elements)[0] + points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) + num_of_elt=np.shape(model.mesh.elements)[0] - num_of_points=numpy.size(points)[0] - dim=numpy.size(points)[1] - point_per_elt=numpy.shape(model.mesh.elements)[1] + num_of_points=np.size(points)[0] + dim=np.size(points)[1] + point_per_elt=np.shape(model.mesh.elements)[1] celltype=5 #triangles @@ -67,6 +67,6 @@ for solution in solnames: #looking for multiple time steps - if (numpy.size(res_struct.__dict__[solution])>num_of_timesteps): - num_of_timesteps=numpy.size(res_struct.__dict__[solution]) + if (np.size(res_struct.__dict__[solution])>num_of_timesteps): + num_of_timesteps=np.size(res_struct.__dict__[solution]) num_of_timesteps=int(num_of_timesteps)+1 else: @@ -106,11 +106,11 @@ for sol in solnames: #dealing with results on different timesteps - if(numpy.size(res_struct.__dict__[sol])>timestep): + if(np.size(res_struct.__dict__[sol])>timestep): timestep = step else: - timestep = numpy.size(res_struct.__dict__[sol]) + timestep = np.size(res_struct.__dict__[sol]) #getting the fields in the solution - if(numpy.size(res_struct.__dict__[sol])>1): + if(np.size(res_struct.__dict__[sol])>1): fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep-1).__dict__) else: @@ -118,15 +118,15 @@ #check which field is a real result and print for field in fieldnames: - if(numpy.size(res_struct.__dict__[sol])>1): + if(np.size(res_struct.__dict__[sol])>1): fieldstruct=res_struct.__dict__[sol].__getitem__(timestep-1).__dict__[field] else: fieldstruct=res_struct.__dict__[sol].__dict__[field] - if ((numpy.size(fieldstruct))==num_of_points): + if ((np.size(fieldstruct))==num_of_points): fid.write('SCALARS %s float 1 \n' % field) fid.write('LOOKUP_TABLE default\n') for node in range(0,num_of_points): #paraview does not like NaN, replacing - if numpy.isnan(fieldstruct[node]): + if np.isnan(fieldstruct[node]): fid.write('%e\n' % -9999.9999) #also checking for verry small value that mess up @@ -143,10 +143,10 @@ for field in othernames: # fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe)); - if ((numpy.size(other_struct.__dict__[field]))==num_of_points): + if ((np.size(other_struct.__dict__[field]))==num_of_points): fid.write('SCALARS %s float 1 \n' % field) fid.write('LOOKUP_TABLE default\n') for node in range(0,num_of_points): #paraview does not like NaN, replacing - if numpy.isnan(other_struct.__dict__[field][node]): + if np.isnan(other_struct.__dict__[field][node]): fid.write('%e\n' % -9999.9999) #also checking for verry small value that mess up Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 21302) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os import model @@ -42,13 +42,13 @@ #get the element related variables if 'z' in dict.keys(model.mesh.__dict__): - points=numpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) + points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) dim=3 else: - points=numpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x)))) + points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) dim=2 - num_of_points=numpy.size(model.mesh.x) - num_of_elt=numpy.shape(model.mesh.elements)[0] - point_per_elt=numpy.shape(model.mesh.elements)[1] + num_of_points=np.size(model.mesh.x) + num_of_elt=np.shape(model.mesh.elements)[0] + point_per_elt=np.shape(model.mesh.elements)[1] #Select the type of element function of the number of nodes per elements @@ -70,6 +70,6 @@ for solution in solnames: #looking for multiple time steps - if (numpy.size(res_struct.__dict__[solution])>num_of_timesteps): - num_of_timesteps=numpy.size(res_struct.__dict__[solution]) + if (np.size(res_struct.__dict__[solution])>num_of_timesteps): + num_of_timesteps=np.size(res_struct.__dict__[solution]) num_of_timesteps=int(num_of_timesteps) else: @@ -112,11 +112,11 @@ for sol in solnames: #dealing with results on different timesteps - if(numpy.size(res_struct.__dict__[sol])>timestep): + if(np.size(res_struct.__dict__[sol])>timestep): timestep = step else: - timestep = numpy.size(res_struct.__dict__[sol]) + timestep = np.size(res_struct.__dict__[sol]) #getting the fields in the solution - if(numpy.size(res_struct.__dict__[sol])>1): + if(np.size(res_struct.__dict__[sol])>1): fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep).__dict__) else: @@ -124,15 +124,15 @@ #check which field is a real result and print for field in fieldnames: - if(numpy.size(res_struct.__dict__[sol])>1): + if(np.size(res_struct.__dict__[sol])>1): fieldstruct=res_struct.__dict__[sol].__getitem__(timestep).__dict__[field] else: fieldstruct=res_struct.__dict__[sol].__dict__[field] - if ((numpy.size(fieldstruct))==num_of_points): + if ((np.size(fieldstruct))==num_of_points): fid.write('SCALARS %s float 1 \n' % field) fid.write('LOOKUP_TABLE default\n') for node in range(0,num_of_points): #paraview does not like NaN, replacing - if numpy.isnan(fieldstruct[node]): + if np.isnan(fieldstruct[node]): fid.write('%e\n' % -9999.9999) #also checking for verry small value that mess up @@ -148,10 +148,10 @@ othernames=(dict.keys(other_struct.__dict__)) for field in othernames: - if ((numpy.size(other_struct.__dict__[field]))==num_of_points): + if ((np.size(other_struct.__dict__[field]))==num_of_points): fid.write('SCALARS %s float 1 \n' % field) fid.write('LOOKUP_TABLE default\n') for node in range(0,num_of_points): #paraview does not like NaN, replacing - if numpy.isnan(other_struct.__dict__[field][node]): + if np.isnan(other_struct.__dict__[field][node]): fid.write('%e\n' % -9999.9999) #also checking for verry small value that mess up Index: /issm/trunk-jpl/src/m/contrib/morlighem/bamg/YamsCall.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/bamg/YamsCall.py (revision 21302) +++ /issm/trunk-jpl/src/m/contrib/morlighem/bamg/YamsCall.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import time import subprocess @@ -37,5 +37,5 @@ t1=time.time() print "%s" % ' computing metric...' - metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,numpy.empty(0,int)) + metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,np.empty(0,int)) t2=time.time() print "%s%d%s\n" % (' done (',t2-t1,' seconds)') @@ -44,5 +44,5 @@ t1=time.time() print "%s" % ' writing initial mesh files...' - numpy.savetxt('carre0.met',metric) + np.savetxt('carre0.met',metric) f=open('carre0.mesh','w') @@ -66,12 +66,12 @@ #Deal with rifts - if numpy.any(not numpy.isnan(md.rifts.riftstruct)): + if np.any(not np.isnan(md.rifts.riftstruct)): #we have the list of triangles that make up the rift. keep those triangles around during refinement. - triangles=numpy.empty(0,int) + triangles=np.empty(0,int) for riftstruct in md.rifts.riftstruct: - triangles=numpy.concatenate((triangles,riftstruct.segments[:,2])) + triangles=np.concatenate((triangles,riftstruct.segments[:,2])) - f.write("\n\n%s\n%i\n\n" % ('RequiredTriangles',numpy.size(triangles))) + f.write("\n\n%s\n%i\n\n" % ('RequiredTriangles',np.size(triangles))) for triangle in triangles: f.write("%i\n" % triangle) @@ -97,12 +97,12 @@ t1=time.time() print "\n%s" % ' reading final mesh files...' - Tria=numpy.loadtxt('carre1.tria',int) - Coor=numpy.loadtxt('carre1.coor',float) + Tria=np.loadtxt('carre1.tria',int) + Coor=np.loadtxt('carre1.coor',float) md.mesh.x=Coor[:,0] md.mesh.y=Coor[:,1] - md.mesh.z=numpy.zeros((numpy.size(Coor,axis=0),1)) + md.mesh.z=np.zeros((np.size(Coor,axis=0),1)) md.mesh.elements=Tria - md.mesh.numberofvertices=numpy.size(Coor,axis=0) - md.mesh.numberofelements=numpy.size(Tria,axis=0) + md.mesh.numberofvertices=np.size(Coor,axis=0) + md.mesh.numberofelements=np.size(Tria,axis=0) numberofelements2=md.mesh.numberofelements t2=time.time() Index: /issm/trunk-jpl/src/m/coordsystems/gmtmask.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gmtmask.py (revision 21302) +++ /issm/trunk-jpl/src/m/coordsystems/gmtmask.py (revision 21303) @@ -1,5 +1,5 @@ from MatlabFuncs import * from model import * -from numpy import * +from np.import * from os import getenv, putenv import subprocess Index: /issm/trunk-jpl/src/m/coordsystems/ll2xy.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/ll2xy.py (revision 21302) +++ /issm/trunk-jpl/src/m/coordsystems/ll2xy.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71): Index: /issm/trunk-jpl/src/m/coordsystems/xy2ll.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/xy2ll.py (revision 21302) +++ /issm/trunk-jpl/src/m/coordsystems/xy2ll.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from math import pi @@ -37,8 +37,8 @@ raise StandardError('bad usage: type "help(xy2ll)" for details') - # if x,y passed as lists, convert to numpy arrays - if type(x) != "numpy.ndarray": + # if x,y passed as lists, convert to np.arrays + if type(x) != "np.ndarray": x=np.array(x) - if type(y) != "numpy.ndarray": + if type(y) != "np.ndarray": y=np.array(y) Index: /issm/trunk-jpl/src/m/exp/expcoarsen.py =================================================================== --- /issm/trunk-jpl/src/m/exp/expcoarsen.py (revision 21302) +++ /issm/trunk-jpl/src/m/exp/expcoarsen.py (revision 21303) @@ -1,4 +1,4 @@ import os.path -import numpy as np +import numpy as np from collections import OrderedDict from expread import expread Index: /issm/trunk-jpl/src/m/exp/expdisp.py =================================================================== --- /issm/trunk-jpl/src/m/exp/expdisp.py (revision 21302) +++ /issm/trunk-jpl/src/m/exp/expdisp.py (revision 21303) @@ -1,4 +1,4 @@ from expread import expread -import numpy as np +import numpy as np from matplotlib.path import Path import matplotlib.patches as patches Index: /issm/trunk-jpl/src/m/exp/expread.py =================================================================== --- /issm/trunk-jpl/src/m/exp/expread.py (revision 21302) +++ /issm/trunk-jpl/src/m/exp/expread.py (revision 21303) @@ -1,9 +1,11 @@ import os.path -import numpy +import numpy as np from collections import OrderedDict import MatlabFuncs as m def expread(filename): + """ + EXPREAD - read a file exp and build a Structure @@ -23,6 +25,6 @@ See also EXPDOC, EXPWRITEASVERTICES + """ - #some checks if not os.path.exists(filename): @@ -31,14 +33,10 @@ #initialize number of profile contours=[] - #open file fid=open(filename,'r') - #loop over the number of profiles while True: - #update number of profiles contour=OrderedDict() - #Get file name A=fid.readline() @@ -50,4 +48,5 @@ if not (len(A) == 2 and m.strcmp(A[0],'##') and m.strncmp(A[1],'Name:',5)): break + if len(A[1])>5: contour['name']=A[1][5:-1] @@ -59,5 +58,4 @@ if not (len(A) == 2 and m.strcmp(A[0],'##') and m.strncmp(A[1],'Icon:',5)): break - #Get Info A=fid.readline().split() @@ -72,11 +70,10 @@ #Get Info A=fid.readline().split() - if not (len(A) == 5 and m.strcmp(A[0],'#') and m.strcmp(A[1],'X') and m.strcmp(A[2],'pos') \ - and m.strcmp(A[3],'Y') and m.strcmp(A[4],'pos')): + if not (len(A) == 5 and m.strcmp(A[0],'#') and m.strcmp(A[1],'X') and m.strcmp(A[2],'pos') + and m.strcmp(A[3],'Y') and m.strcmp(A[4],'pos')): break - #Get Coordinates - contour['x']=numpy.empty(contour['nods']) - contour['y']=numpy.empty(contour['nods']) + contour['x']=np.empty(contour['nods']) + contour['y']=np.empty(contour['nods']) for i in xrange(int(contour['nods'])): A=fid.readline().split() @@ -93,8 +90,5 @@ contours.append(contour) - #close file fid.close() - return contours - Index: /issm/trunk-jpl/src/m/exp/expwrite.py =================================================================== --- /issm/trunk-jpl/src/m/exp/expwrite.py (revision 21302) +++ /issm/trunk-jpl/src/m/exp/expwrite.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def expwrite(contours,filename): @@ -21,5 +21,5 @@ fid=open(filename,'w') for x,y in zip(contours['x'],contours['y']): - #if numpy.size(contour['x'])!=numpy.size(contour['y']): + #if np.size(contour['x'])!=np.size(contour['y']): if len(x)!=len(y): raise RuntimeError("contours x and y coordinates must be of identical size") @@ -36,6 +36,6 @@ fid.write("%s\n" % '## Icon:0') fid.write("%s\n" % '# Points Count Value') - #fid.write("%i %f\n" % (numpy.size(contour['x']),contour['density'])) - fid.write("%i %f\n" % (numpy.size(x),density)) + #fid.write("%i %f\n" % (np.size(contour['x']),contour['density'])) + fid.write("%i %f\n" % (np.size(x),density)) fid.write("%s\n" % '# X pos Y pos') #for x,y in zip(contour['x'],contour['y']): Index: /issm/trunk-jpl/src/m/extrusion/DepthAverage.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/DepthAverage.py (revision 21302) +++ /issm/trunk-jpl/src/m/extrusion/DepthAverage.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from project2d import project2d @@ -46,5 +46,5 @@ if vec2d: - vector_average=vector_average.reshape(-1,1) + vector_average=vector_average.reshape(-1,) return vector_average Index: /issm/trunk-jpl/src/m/extrusion/project2d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project2d.py (revision 21302) +++ /issm/trunk-jpl/src/m/extrusion/project2d.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np def project2d(md3d,value,layer): @@ -45,5 +45,5 @@ if vec2d: - projection_value=projection_value.reshape(-1,1) + projection_value=projection_value.reshape(-1,) return projection_value Index: /issm/trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 21302) +++ /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from pairoptions import pairoptions @@ -37,9 +37,9 @@ vector1d=False - if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1: + if isinstance(vector2d,np.ndarray) and np.ndim(vector2d)==1: vector1d=True - vector2d=vector2d.reshape(-1,1) + vector2d=vector2d.reshape(-1,) - if isinstance(vector2d,(bool,int,long,float)) or numpy.size(vector2d)==1: + if isinstance(vector2d,(bool,int,long,float)) or np.size(vector2d)==1: projected_vector=vector2d @@ -47,38 +47,69 @@ #Initialize 3d vector - if vector2d.shape[0]==md.mesh.numberofvertices2d: - projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: - projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - projected_vector[-1,:]=vector2d[-1,:] - vector2d=vector2d[:-1,:] + if np.ndim(vector2d)==1: + if vector2d.shape[0]==md.mesh.numberofvertices2d: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype) + elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1))).astype(vector2d.dtype) + projected_vector[-1]=vector2d[-1] + vector2d=vector2d[:-1] + else: + raise TypeError("vector length not supported") + #Fill in + if layer==0: + for i in xrange(md.mesh.numberoflayers): + projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d)]=vector2d + else: + projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d)]=vector2d else: - raise TypeError("vector length not supported") + if vector2d.shape[0]==md.mesh.numberofvertices2d: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices,np.size(vector2d,axis=1)))).astype(vector2d.dtype) + elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) + projected_vector[-1,:]=vector2d[-1,:] + vector2d=vector2d[:-1,:] + else: + raise TypeError("vector length not supported") + #Fill in + if layer==0: + for i in xrange(md.mesh.numberoflayers): + projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d),:]=vector2d + else: + projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d - #Fill in - if layer==0: - for i in xrange(md.mesh.numberoflayers): - projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d),:]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d elif vectype.lower()=='element': #Initialize 3d vector - if vector2d.shape[0]==md.mesh.numberofelements2d: - projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif vector2d.shape[0]==md.mesh.numberofelements2d+1: - projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - projected_vector[-1,:]=vector2d[-1,:] - vector2d=vector2d[:-1,:] + if np.ndim(vector2d)==1: + if vector2d.shape[0]==md.mesh.numberofelements2d: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements))).astype(vector2d.dtype) + elif vector2d.shape[0]==md.mesh.numberofelements2d+1: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1))).astype(vector2d.dtype) + projected_vector[-1]=vector2d[-1] + vector2d=vector2d[:-1] + else: + raise TypeError("vector length not supported") + #Fill in + if layer==0: + for i in xrange(md.mesh.numberoflayers-1): + projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d)]=vector2d + else: + projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d)]=vector2d else: - raise TypeError("vector length not supported") - - #Fill in - if layer==0: - for i in xrange(md.mesh.numberoflayers-1): - projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d),:]=vector2d - else: - projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d),:]=vector2d + if vector2d.shape[0]==md.mesh.numberofelements2d: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements, np.size(vector2d,axis=1)))).astype(vector2d.dtype) + elif vector2d.shape[0]==md.mesh.numberofelements2d+1: + projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) + projected_vector[-1,:]=vector2d[-1,:] + vector2d=vector2d[:-1,:] + else: + raise TypeError("vector length not supported") + #Fill in + if layer==0: + for i in xrange(md.mesh.numberoflayers-1): + projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d),:]=vector2d + else: + projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d),:]=vector2d else: Index: /issm/trunk-jpl/src/m/geometry/FlagElements.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/FlagElements.py (revision 21302) +++ /issm/trunk-jpl/src/m/geometry/FlagElements.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os #from basinzoom import basinzoon @@ -24,8 +24,8 @@ if isinstance(region,(str,unicode)): if not region: - flag=numpy.zeros(md.mesh.numberofelements,bool) + flag=np.zeros(md.mesh.numberofelements,bool) invert=0 elif m.strcmpi(region,'all'): - flag=numpy.ones(md.mesh.numberofelements,bool) + flag=np.ones(md.mesh.numberofelements,bool) invert=0 else: @@ -44,5 +44,5 @@ xlim,ylim=basinzoom('basin',region) flag_nodes=p.logical_and_n(md.mesh.xxlim[0],md.mesh.yylim[0]) - flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool) + flag=np.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool) else: #ok, flag elements @@ -51,11 +51,11 @@ if invert: - flag=numpy.logical_not(flag) + flag=np.logical_not(flag) - elif isinstance(region,numpy.ndarray) or isinstance(region,bool): - if numpy.size(region,0)==md.mesh.numberofelements: + elif isinstance(region,np.ndarray) or isinstance(region,bool): + if np.size(region,0)==md.mesh.numberofelements: flag=region - elif numpy.size(region,0)==md.mesh.numberofvertices: - flag=(numpy.sum(region[md.mesh.elements-1]>0,axis=1)==numpy.size(md.mesh.elements,1)) + elif np.size(region,0)==md.mesh.numberofvertices: + flag=(np.sum(region[md.mesh.elements-1]>0,axis=1)==np.size(md.mesh.elements,1)) else: raise TypeError("Flaglist for region must be of same size as number of elements in model.") Index: /issm/trunk-jpl/src/m/geometry/GetAreas.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/GetAreas.py (revision 21302) +++ /issm/trunk-jpl/src/m/geometry/GetAreas.py (revision 21303) @@ -1,5 +1,5 @@ -import numpy +import numpy as np -def GetAreas(index,x,y,z=numpy.array([])): +def GetAreas(index,x,y,z=np.array([])): """ GETAREAS - compute areas or volumes of elements @@ -18,19 +18,19 @@ #get number of elements and number of nodes - nels=numpy.size(index,axis=0) - nods=numpy.size(x) + nels=np.size(index,axis=0) + nods=np.size(x) #some checks - if numpy.size(y)!=nods or (z and numpy.size(z)!=nods): + if np.size(y)!=nods or (z and np.size(z)!=nods): raise TypeError("GetAreas error message: x,y and z do not have the same length.") - if numpy.max(index)>nods: + if np.max(index)>nods: raise TypeError("GetAreas error message: index should not have values above %d." % nods) - if (not z and numpy.size(index,axis=1)!=3): + if (not z and np.size(index,axis=1)!=3): raise TypeError("GetAreas error message: index should have 3 columns for 2d meshes.") - if (z and numpy.size(index,axis=1)!=6): + if (z and np.size(index,axis=1)!=6): raise TypeError("GetAreas error message: index should have 6 columns for 3d meshes.") #initialization - areas=numpy.zeros(nels) + areas=np.zeros(nels) x1=x[index[:,0]-1] x2=x[index[:,1]-1] @@ -46,5 +46,5 @@ else: #V=area(triangle)*1/3(z1+z2+z3) - thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1]) + thickness=np.mean(z[index[:,3:6]-1])-np.mean(z[index[:,0:3]-1]) areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness Index: /issm/trunk-jpl/src/m/geometry/SegIntersect.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/SegIntersect.py (revision 21302) +++ /issm/trunk-jpl/src/m/geometry/SegIntersect.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def SegIntersect(seg1,seg2): @@ -24,14 +24,14 @@ yD=seg2[1,1] - O2A=numpy.array([xA,yA])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.]) - O2B=numpy.array([xB,yB])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.]) - O1C=numpy.array([xC,yC])-numpy.array([xA/2.+xB/2.,yB/2.+yA/2.]) - O1D=numpy.array([xD,yD])-numpy.array([xA/2.+xB/2.,yB/2.+yA/2.]) + O2A=np.array([xA,yA])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) + O2B=np.array([xB,yB])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) + O1C=np.array([xC,yC])-np.array([xA/2.+xB/2.,yB/2.+yA/2.]) + O1D=np.array([xD,yD])-np.array([xA/2.+xB/2.,yB/2.+yA/2.]) - n1=numpy.array([yA-yB,xB-xA]) #normal vector to segA - n2=numpy.array([yC-yD,xD-xC]) #normal vector to segB + n1=np.array([yA-yB,xB-xA]) #normal vector to segA + n2=np.array([yC-yD,xD-xC]) #normal vector to segB - test1=numpy.dot(n2,O2A) - test2=numpy.dot(n2,O2B) + test1=np.dot(n2,O2A) + test2=np.dot(n2,O2B) if test1*test2>0: @@ -39,6 +39,6 @@ return bval - test3=numpy.dot(n1,O1C) - test4=numpy.dot(n1,O1D) + test3=np.dot(n1,O1C) + test4=np.dot(n1,O1D) if test3*test4>0: @@ -47,12 +47,12 @@ #if colinear - if test1*test2==0 and test3*test4==0 and numpy.linalg.det(numpy.hstack((n1.reshape((-1,1)),n2.reshape(-1,1))))==0: + if test1*test2==0 and test3*test4==0 and np.linalg.det(np.hstack((n1.reshape((-1,)),n2.reshape(-1,))))==0: #projection on the axis O1O2 - O2O1=numpy.array([xA/2.+xB/2.,yB/2.+yA/2.])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.]) - O1A=numpy.dot(O2O1,(O2A-O2O1)) - O1B=numpy.dot(O2O1,(O2B-O2O1)) - O1C=numpy.dot(O2O1,O1C) - O1D=numpy.dot(O2O1,O1D) + O2O1=np.array([xA/2.+xB/2.,yB/2.+yA/2.])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) + O1A=np.dot(O2O1,(O2A-O2O1)) + O1B=np.dot(O2O1,(O2B-O2O1)) + O1C=np.dot(O2O1,O1C) + O1D=np.dot(O2O1,O1D) #test if one point is included in the other segment (->bval=1) Index: /issm/trunk-jpl/src/m/geometry/slope.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/slope.py (revision 21302) +++ /issm/trunk-jpl/src/m/geometry/slope.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff Index: /issm/trunk-jpl/src/m/interp/SectionValues.py =================================================================== --- /issm/trunk-jpl/src/m/interp/SectionValues.py (revision 21302) +++ /issm/trunk-jpl/src/m/interp/SectionValues.py (revision 21303) @@ -1,5 +1,5 @@ import os from expread import expread -import numpy as np +import numpy as np from project2d import project2d #from InterpFromMesh2d import InterpFromMesh2d Index: /issm/trunk-jpl/src/m/interp/averaging.py =================================================================== --- /issm/trunk-jpl/src/m/interp/averaging.py (revision 21302) +++ /issm/trunk-jpl/src/m/interp/averaging.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from GetAreas import GetAreas import MatlabFuncs as m @@ -82,5 +82,5 @@ average_node = csc_matrix(average_node) else: - average_node=csc_matrix(data.reshape(-1,1)) + average_node=csc_matrix(data.reshape(-1,)) #loop over iteration Index: /issm/trunk-jpl/src/m/interp/holefiller.py =================================================================== --- /issm/trunk-jpl/src/m/interp/holefiller.py (revision 21302) +++ /issm/trunk-jpl/src/m/interp/holefiller.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from scipy.spatial import cKDTree Index: /issm/trunk-jpl/src/m/interp/interp.py =================================================================== --- /issm/trunk-jpl/src/m/interp/interp.py (revision 21302) +++ /issm/trunk-jpl/src/m/interp/interp.py (revision 21303) @@ -1,4 +1,4 @@ # module for inperpolating/smoothing data -import numpy as np +import numpy as np from scipy.interpolate import CloughTocher2DInterpolator, Rbf from scipy.spatial import cKDTree Index: /issm/trunk-jpl/src/m/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadvars.py (revision 21302) +++ /issm/trunk-jpl/src/m/io/loadvars.py (revision 21303) @@ -1,5 +1,5 @@ import shelve import os.path -import numpy as np +import numpy as np from netCDF4 import Dataset from netCDF4 import chartostring @@ -7,4 +7,5 @@ from os import path from whichdb import whichdb +from collections import OrderedDict from model import * @@ -130,5 +131,5 @@ Tree[t].__dict__[str(var)]=varval[0] else: - if varval[0]==u'': + if str(varval[0])=='': Tree.__dict__[str(var)]=[] elif varval[0]=='True': @@ -167,5 +168,5 @@ #dealling with dict if varval.dtype==str: - Tree.__dict__[str(var)]=dict(zip(varval[:,0], varval[:,1])) + Tree.__dict__[str(var)]=OrderedDict(zip(varval[:,0], varval[:,1])) else: if type(Tree)==list: @@ -178,4 +179,8 @@ else: Tree.__dict__[str(var)]=varval[:,:] + try: + print('treated {}.{} with type {}'.format(mod,var,type(Tree.__dict__[str(var)]))) + except AttributeError: + print('treated {}.{} with type {}'.format(mod,var,type(Tree[0].__dict__[str(var)]))) elif vardim==3: if type(Tree)==list: Index: /issm/trunk-jpl/src/m/materials/TMeltingPoint.py =================================================================== --- /issm/trunk-jpl/src/m/materials/TMeltingPoint.py (revision 21302) +++ /issm/trunk-jpl/src/m/materials/TMeltingPoint.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np def TMeltingPoint(reftemp,pressure): Index: /issm/trunk-jpl/src/m/materials/cuffey.py =================================================================== --- /issm/trunk-jpl/src/m/materials/cuffey.py (revision 21302) +++ /issm/trunk-jpl/src/m/materials/cuffey.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def cuffey(temperature): @@ -13,35 +13,35 @@ """ - if numpy.any(temperature<0.): + if np.any(temperature<0.): raise RuntimeError("input temperature should be in Kelvin (positive)") T = temperature-273.15 - rigidity=numpy.zeros_like(T) - pos=numpy.nonzero(T<=-45) + rigidity=np.zeros_like(T) + pos=np.nonzero(T<=-45) rigidity[pos]=10**8*(-0.000396645116301*(T[pos]+50)**3+ 0.013345579471334*(T[pos]+50)**2 -0.356868703259105*(T[pos]+50)+7.272363035371383) - pos=numpy.nonzero(numpy.logical_and(-45<=T,T<-40)) + pos=np.nonzero(np.logical_and(-45<=T,T<-40)) rigidity[pos]=10**8*(-0.000396645116301*(T[pos]+45)**3+ 0.007395902726819*(T[pos]+45)**2 -0.253161292268336*(T[pos]+45)+5.772078366321591) - pos=numpy.nonzero(numpy.logical_and(-40<=T,T<-35)) + pos=np.nonzero(np.logical_and(-40<=T,T<-35)) rigidity[pos]=10**8*(0.000408322072669*(T[pos]+40)**3+ 0.001446225982305*(T[pos]+40)**2 -0.208950648722716*(T[pos]+40)+4.641588833612773) - pos=numpy.nonzero(numpy.logical_and(-35<=T,T<-30)) + pos=np.nonzero(np.logical_and(-35<=T,T<-30)) rigidity[pos]=10**8*(-0.000423888728124*(T[pos]+35)**3+ 0.007571057072334*(T[pos]+35)**2 -0.163864233449525*(T[pos]+35)+3.684031498640382) - pos=numpy.nonzero(numpy.logical_and(-30<=T,T<-25)) + pos=np.nonzero(np.logical_and(-30<=T,T<-25)) rigidity[pos]=10**8*(0.000147154327025*(T[pos]+30)**3+ 0.001212726150476*(T[pos]+30)**2 -0.119945317335478*(T[pos]+30)+3.001000667185614) - pos=numpy.nonzero(numpy.logical_and(-25<=T,T<-20)) + pos=np.nonzero(np.logical_and(-25<=T,T<-20)) rigidity[pos]=10**8*(-0.000193435838672*(T[pos]+25)**3+ 0.003420041055847*(T[pos]+25)**2 -0.096781481303861*(T[pos]+25)+2.449986525148220) - pos=numpy.nonzero(numpy.logical_and(-20<=T,T<-15)) + pos=np.nonzero(np.logical_and(-20<=T,T<-15)) rigidity[pos]=10**8*(0.000219771255067*(T[pos]+20)**3+ 0.000518503475772*(T[pos]+20)**2 -0.077088758645767*(T[pos]+20)+2.027400665191131) - pos=numpy.nonzero(numpy.logical_and(-15<=T,T<-10)) + pos=np.nonzero(np.logical_and(-15<=T,T<-10)) rigidity[pos]=10**8*(-0.000653438900191*(T[pos]+15)**3+ 0.003815072301777*(T[pos]+15)**2 -0.055420879758021*(T[pos]+15)+1.682390865739973) - pos=numpy.nonzero(numpy.logical_and(-10<=T,T<-5)) + pos=np.nonzero(np.logical_and(-10<=T,T<-5)) rigidity[pos]=10**8*(0.000692439419762*(T[pos]+10)**3 -0.005986511201093 *(T[pos]+10)**2 -0.066278074254598*(T[pos]+10)+1.418983411970382) - pos=numpy.nonzero(numpy.logical_and(-5<=T,T<-2)) + pos=np.nonzero(np.logical_and(-5<=T,T<-2)) rigidity[pos]=10**8*(-0.000132282004110*(T[pos]+5)**3 +0.004400080095332*(T[pos]+5)**2 -0.074210229783403*(T[pos]+5)+ 1.024485188140279) - pos=numpy.nonzero(-2<=T) + pos=np.nonzero(-2<=T) rigidity[pos]=10**8*(-0.000132282004110*(T[pos]+2)**3 +0.003209542058346*(T[pos]+2)**2 -0.051381363322371*(T[pos]+2)+ 0.837883605537096) #Now make sure that rigidity is positive - pos=numpy.nonzero(rigidity<0) + pos=np.nonzero(rigidity<0) rigidity[pos]=1**6 Index: /issm/trunk-jpl/src/m/materials/cuffeytemperate.py =================================================================== --- /issm/trunk-jpl/src/m/materials/cuffeytemperate.py (revision 21302) +++ /issm/trunk-jpl/src/m/materials/cuffeytemperate.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import cuffey @@ -15,15 +15,15 @@ """ - if numpy.any(temperature<0.): + if np.any(temperature<0.): raise RuntimeError("input temperature should be in Kelvin (positive)") - if (numpy.any(temperature.shape~=waterfraction.shape)), + if (np.any(temperature.shape~=waterfraction.shape)), error('input temperature and waterfraction should have same size!'); end - if numpy.any(waterfraction<0 | waterfraction>1) + if np.any(waterfraction<0 | waterfraction>1) error('input waterfraction should be between 0 and 1'); end - rigidity=numpy.multiply(cuffey(temperature), (1*numpy.ones(waterfraction.shape)+181.25*numpy.maximum(numpy.zeros(waterfraction.shape), numpy.minimum(0.01*numpy.ones(waterfraction.shape), waterfraction)))**(-1/stressexp)); + rigidity=np.multiply(cuffey(temperature), (1*np.ones(waterfraction.shape)+181.25*np.maximum(np.zeros(waterfraction.shape), np.minimum(0.01*np.ones(waterfraction.shape), waterfraction)))**(-1/stressexp)); return rigidity Index: /issm/trunk-jpl/src/m/materials/paterson.py =================================================================== --- /issm/trunk-jpl/src/m/materials/paterson.py (revision 21302) +++ /issm/trunk-jpl/src/m/materials/paterson.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def paterson(temperature): @@ -12,12 +12,12 @@ """ - if numpy.any(temperature<0.): + if np.any(temperature<0.): raise RuntimeError("input temperature should be in Kelvin (positive)") - if numpy.ndim(temperature)==2: + if np.ndim(temperature)==2: #T = temperature.reshape(-1,)-273.15 T = temperature.flatten()-273.15 elif isinstance(temperature,float) or isinstance(temperature,int): - T = numpy.array([temperature])-273.15 + T = np.array([temperature])-273.15 else: T = temperature-273.15 @@ -36,41 +36,41 @@ # rigidity=fittedmodel(temperature); - rigidity=numpy.zeros_like(T) - pos1=numpy.nonzero(T<=-45)[0] + rigidity=np.zeros_like(T) + pos1=np.nonzero(T<=-45)[0] if len(pos1): rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) - pos2=numpy.nonzero(numpy.logical_and(-45<=T,T<-40))[0] + pos2=np.nonzero(np.logical_and(-45<=T,T<-40))[0] if len(pos2): rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) - pos3=numpy.nonzero(numpy.logical_and(-40<=T,T<-35))[0] + pos3=np.nonzero(np.logical_and(-40<=T,T<-35))[0] if len(pos3): rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) - pos4=numpy.nonzero(numpy.logical_and(-35<=T,T<-30))[0] + pos4=np.nonzero(np.logical_and(-35<=T,T<-30))[0] if len(pos4): rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) - pos5=numpy.nonzero(numpy.logical_and(-30<=T,T<-25))[0] + pos5=np.nonzero(np.logical_and(-30<=T,T<-25))[0] if len(pos5): rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) - pos6=numpy.nonzero(numpy.logical_and(-25<=T,T<-20))[0] + pos6=np.nonzero(np.logical_and(-25<=T,T<-20))[0] if len(pos6): rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) - pos7=numpy.nonzero(numpy.logical_and(-20<=T,T<-15))[0] + pos7=np.nonzero(np.logical_and(-20<=T,T<-15))[0] if len(pos7): rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) - pos8=numpy.nonzero(numpy.logical_and(-15<=T,T<-10))[0] + pos8=np.nonzero(np.logical_and(-15<=T,T<-10))[0] if len(pos8): rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) - pos9=numpy.nonzero(numpy.logical_and(-10<=T,T<-5))[0] + pos9=np.nonzero(np.logical_and(-10<=T,T<-5))[0] if len(pos9): rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) - pos10=numpy.nonzero(numpy.logical_and(-5<=T,T<-2))[0] + pos10=np.nonzero(np.logical_and(-5<=T,T<-2))[0] if len(pos10): rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) - pos11=numpy.nonzero(-2<=T)[0] + pos11=np.nonzero(-2<=T)[0] if len(pos11): rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) #Now make sure that rigidity is positive - pos=numpy.nonzero(rigidity<0)[0] + pos=np.nonzero(rigidity<0)[0] if len(pos): rigidity[pos]=1.e6 Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.py =================================================================== --- /issm/trunk-jpl/src/m/mech/analyticaldamage.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/analyticaldamage.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from averaging import averaging #from plotmodel import plotmodel Index: /issm/trunk-jpl/src/m/mech/backstressfrominversion.py =================================================================== --- /issm/trunk-jpl/src/m/mech/backstressfrominversion.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/backstressfrominversion.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from averaging import averaging from thomasparams import thomasparams Index: /issm/trunk-jpl/src/m/mech/calcbackstress.py =================================================================== --- /issm/trunk-jpl/src/m/mech/calcbackstress.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/calcbackstress.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from averaging import averaging from thomasparams import thomasparams Index: /issm/trunk-jpl/src/m/mech/damagefrominversion.py =================================================================== --- /issm/trunk-jpl/src/m/mech/damagefrominversion.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/damagefrominversion.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np def damagefrominversion(md): Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.py =================================================================== --- /issm/trunk-jpl/src/m/mech/mechanicalproperties.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff from results import results Index: /issm/trunk-jpl/src/m/mech/robintemperature.py =================================================================== --- /issm/trunk-jpl/src/m/mech/robintemperature.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/robintemperature.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from scipy.special import erf Index: /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py =================================================================== --- /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np def steadystateiceshelftemp(md,surfacetemp,basaltemp): @@ -18,5 +18,4 @@ temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp) """ - if len(md.geometry.thickness)!=md.mesh.numberofvertices: raise ValueError('steadystateiceshelftemp error message: thickness should have a length of ' + md.mesh.numberofvertices) Index: /issm/trunk-jpl/src/m/mech/thomasparams.py =================================================================== --- /issm/trunk-jpl/src/m/mech/thomasparams.py (revision 21302) +++ /issm/trunk-jpl/src/m/mech/thomasparams.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy as np +import numpy as np from averaging import averaging Index: /issm/trunk-jpl/src/m/mesh/ComputeHessian.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/ComputeHessian.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/ComputeHessian.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff from GetAreas import GetAreas @@ -20,9 +20,9 @@ #some variables - numberofnodes=numpy.size(x) - numberofelements=numpy.size(index,axis=0) + numberofnodes=np.size(x) + numberofelements=np.size(index,axis=0) #some checks - if numpy.size(field)!=numberofnodes and numpy.size(field)!=numberofelements: + if np.size(field)!=numberofnodes and np.size(field)!=numberofelements: raise TypeError("ComputeHessian error message: the given field size not supported yet") if not m.strcmpi(type,'node') and not m.strcmpi(type,'element'): @@ -38,28 +38,28 @@ #compute weights that hold the volume of all the element holding the node i - weights=m.sparse(line,numpy.ones((linesize,1)),numpy.tile(areas.reshape(-1,1),(3,1)),numberofnodes,1) + weights=m.sparse(line,np.ones((linesize,1)),np.tile(areas.reshape(-1,),(3,1)),numberofnodes,1) #compute field on nodes if on elements - if numpy.size(field,axis=0)==numberofelements: - field=m.sparse(line,numpy.ones((linesize,1)),numpy.tile(areas*field,(3,1)),numberofnodes,1)/weights + if np.size(field,axis=0)==numberofelements: + field=m.sparse(line,np.ones((linesize,1)),np.tile(areas*field,(3,1)),numberofnodes,1)/weights #Compute gradient for each element - grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1) - grad_ely=numpy.sum(field[index-1,0]*beta,axis=1) + grad_elx=np.sum(field[index-1,0]*alpha,axis=1) + grad_ely=np.sum(field[index-1,0]*beta,axis=1) #Compute gradient for each node (average of the elements around) - gradx=m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*grad_elx).reshape(-1,1),(3,1)),numberofnodes,1) - grady=m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*grad_ely).reshape(-1,1),(3,1)),numberofnodes,1) + gradx=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_elx).reshape(-1,),(3,1)),numberofnodes,1) + grady=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_ely).reshape(-1,),(3,1)),numberofnodes,1) gradx=gradx/weights grady=grady/weights #Compute hessian for each element - hessian=numpy.hstack((numpy.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*beta,axis=1).reshape(-1,1))) + hessian=np.vstack((np.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*beta,axis=1).reshape(-1,))).T if m.strcmpi(type,'node'): #Compute Hessian on the nodes (average of the elements around) - hessian=numpy.hstack((m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,0]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \ - m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,1]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \ - m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,2]).reshape(-1,1),(3,1)),numberofnodes,1)/weights )) + hessian=np.hstack((m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,0]).reshape(-1,),(3,1)),numberofnodes,1)/weights, + m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,1]).reshape(-1,),(3,1)),numberofnodes,1)/weights, + m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,2]).reshape(-1,),(3,1)),numberofnodes,1)/weights )) return hessian Index: /issm/trunk-jpl/src/m/mesh/ComputeMetric.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/ComputeMetric.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/ComputeMetric.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos): @@ -17,21 +17,21 @@ b=hessian[:,1] d=hessian[:,2] - lambda1=0.5*((a+d)+numpy.sqrt(4.*b**2+(a-d)**2)) - lambda2=0.5*((a+d)-numpy.sqrt(4.*b**2+(a-d)**2)) - pos1=numpy.nonzero(lambda1==0.)[0] - pos2=numpy.nonzero(lambda2==0.)[0] - pos3=numpy.nonzero(numpy.logical_and(b==0.,lambda1==lambda2))[0] + lambda1=0.5*((a+d)+np.sqrt(4.*b**2+(a-d)**2)) + lambda2=0.5*((a+d)-np.sqrt(4.*b**2+(a-d)**2)) + pos1=np.nonzero(lambda1==0.)[0] + pos2=np.nonzero(lambda2==0.)[0] + pos3=np.nonzero(np.logical_and(b==0.,lambda1==lambda2))[0] #Modify the eigen values to control the shape of the elements - lambda1=numpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) - lambda2=numpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) + lambda1=np.minimum(np.maximum(np.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) + lambda2=np.minimum(np.maximum(np.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) #compute eigen vectors - norm1=numpy.sqrt(8.*b**2+2.*(d-a)**2+2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2)) + norm1=np.sqrt(8.*b**2+2.*(d-a)**2+2.*(d-a)*np.sqrt((a-d)**2+4.*b**2)) v1x=2.*b/norm1 - v1y=((d-a)+numpy.sqrt((a-d)**2+4.*b**2))/norm1 - norm2=numpy.sqrt(8.*b**2+2.*(d-a)**2-2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2)) + v1y=((d-a)+np.sqrt((a-d)**2+4.*b**2))/norm1 + norm2=np.sqrt(8.*b**2+2.*(d-a)**2-2.*(d-a)*np.sqrt((a-d)**2+4.*b**2)) v2x=2.*b/norm2 - v2y=((d-a)-numpy.sqrt((a-d)**2+4.*b**2))/norm2 + v2y=((d-a)-np.sqrt((a-d)**2+4.*b**2))/norm2 v1x[pos3]=1. @@ -41,32 +41,32 @@ #Compute new metric (for each node M=V*Lambda*V^-1) - metric=numpy.hstack((((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v2y*v1x-lambda2*v1y*v2x)).reshape(-1,1), \ - ((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v1y*v2y-lambda2*v1y*v2y)).reshape(-1,1), \ - ((v1x*v2y-v1y*v2x)**(-1)*(-lambda1*v2x*v1y+lambda2*v1x*v2y)).reshape(-1,1))) + metric=np.vstack((((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v2y*v1x-lambda2*v1y*v2x)).reshape(-1,), + ((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v1y*v2y-lambda2*v1y*v2y)).reshape(-1,), + ((v1x*v2y-v1y*v2x)**(-1)*(-lambda1*v2x*v1y+lambda2*v1x*v2y)).reshape(-1,))).T #some corrections for 0 eigen values - metric[pos1,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos1),1)) - metric[pos2,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos2),1)) + metric[pos1,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos1),1)) + metric[pos2,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos2),1)) #take care of water elements - metric[pos ,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos ),1)) + metric[pos ,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos ),1)) #take care of NaNs if any (use Numpy eig in a loop) - pos=numpy.nonzero(numpy.isnan(metric))[0] - if numpy.size(pos): - print(" %i NaN found in the metric. Use Numpy routine..." % numpy.size(pos)) + pos=np.nonzero(np.isnan(metric))[0] + if np.size(pos): + print(" %i NaN found in the metric. Use Numpy routine..." % np.size(pos)) for posi in pos: - H=numpy.array([[hessian[posi,0],hessian[posi,1]],[hessian[posi,1],hessian[posi,2]]]) - [v,u]=numpy.linalg.eig(H) - v=numpy.diag(v) + H=np.array([[hessian[posi,0],hessian[posi,1]],[hessian[posi,1],hessian[posi,2]]]) + [v,u]=np.linalg.eig(H) + v=np.diag(v) lambda1=v[0,0] lambda2=v[1,1] - v[0,0]=numpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) - v[1,1]=numpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) + v[0,0]=np.minimum(np.maximum(np.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) + v[1,1]=np.minimum(np.maximum(np.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) - metricTria=numpy.dot(numpy.dot(u,v),numpy.linalg.inv(u)) - metric[posi,:]=numpy.array([metricTria[0,0],metricTria[0,1],metricTria[1,1]]) + metricTria=np.dot(np.dot(u,v),np.linalg.inv(u)) + metric[posi,:]=np.array([metricTria[0,0],metricTria[0,1],metricTria[1,1]]) - if numpy.any(numpy.isnan(metric)): + if np.any(np.isnan(metric)): raise RunTimeError("ComputeMetric error message: NaN in the metric despite our efforts...") Index: /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import PythonFuncs as p @@ -12,11 +12,11 @@ """ - edgeelements=numpy.nonzero(\ - p.logical_or_n(numpy.logical_and(elements[:,0]==A,elements[:,1]==B), \ - numpy.logical_and(elements[:,0]==A,elements[:,2]==B), \ - numpy.logical_and(elements[:,1]==A,elements[:,2]==B), \ - numpy.logical_and(elements[:,1]==A,elements[:,0]==B), \ - numpy.logical_and(elements[:,2]==A,elements[:,0]==B), \ - numpy.logical_and(elements[:,2]==A,elements[:,1]==B), \ + edgeelements=np.nonzero(\ + p.logical_or_n(np.logical_and(elements[:,0]==A,elements[:,1]==B), \ + np.logical_and(elements[:,0]==A,elements[:,2]==B), \ + np.logical_and(elements[:,1]==A,elements[:,2]==B), \ + np.logical_and(elements[:,1]==A,elements[:,0]==B), \ + np.logical_and(elements[:,2]==A,elements[:,0]==B), \ + np.logical_and(elements[:,2]==A,elements[:,1]==B), \ ))[0]+1 Index: /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def GetNodalFunctionsCoeff(index,x,y): @@ -23,18 +23,18 @@ #get nels and nods - nels=numpy.size(index,axis=0) - nods=numpy.size(x) + nels=np.size(index,axis=0) + nods=np.size(x) #some checks - if numpy.size(y)!=nods: + if np.size(y)!=nods: raise TypeError("GetNodalFunctionsCoeff error message: x and y do not have the same length.") - if numpy.max(index)>nods: + if np.max(index)>nods: raise TypeError("GetNodalFunctionsCoeff error message: index should not have values above %d." % nods) - if numpy.size(index,axis=1)!=3: + if np.size(index,axis=1)!=3: raise TypeError("GetNodalFunctionsCoeff error message: only 2d meshes supported. index should have 3 columns.") #initialize output - alpha=numpy.zeros((nels,3)) - beta=numpy.zeros((nels,3)) + alpha=np.zeros((nels,3)) + beta=np.zeros((nels,3)) #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma @@ -48,10 +48,10 @@ #get alpha and beta - alpha=numpy.hstack(((invdet*(y2-y3)).reshape(-1,1),(invdet*(y3-y1)).reshape(-1,1),(invdet*(y1-y2)).reshape(-1,1))) - beta =numpy.hstack(((invdet*(x3-x2)).reshape(-1,1),(invdet*(x1-x3)).reshape(-1,1),(invdet*(x2-x1)).reshape(-1,1))) + alpha=np.vstack(((invdet*(y2-y3)).reshape(-1,),(invdet*(y3-y1)).reshape(-1,),(invdet*(y1-y2)).reshape(-1,))).T + beta =np.vstack(((invdet*(x3-x2)).reshape(-1,),(invdet*(x1-x3)).reshape(-1,),(invdet*(x2-x1)).reshape(-1,))).T #get gamma if requested - gamma=numpy.zeros((nels,3)) - gamma=numpy.hstack(((invdet*(x2*y3-x3*y2)).reshape(-1,1),(invdet*(y1*x3-y3*x1)).reshape(-1,1),(invdet*(x1*y2-x2*y1)).reshape(-1,1))) + gamma=np.zeros((nels,3)) + gamma=np.vstack(((invdet*(x2*y3-x3*y2)).reshape(-1,),(invdet*(y1*x3-y3*x1)).reshape(-1,),(invdet*(x1*y2-x2*y1)).reshape(-1,))).T return alpha,beta,gamma Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 21303) @@ -1,4 +1,4 @@ import os.path -import numpy +import numpy as np from mesh2d import mesh2d from collections import OrderedDict @@ -94,13 +94,18 @@ if i: flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0] - if numpy.any(numpy.logical_not(flags)): + if np.any(np.logical_not(flags)): raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") #Add all points to bamg_geometry nods=domaini['nods']-1 #the domain are closed 0=end - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((domaini['x'][0:nods].reshape(-1,1),domaini['y'][0:nods].reshape(-1,1),numpy.ones((nods,1)))))) - bamg_geometry.Edges =numpy.vstack((bamg_geometry.Edges, numpy.hstack((numpy.arange(count+1,count+nods+1).reshape(-1,1),numpy.hstack((numpy.arange(count+2,count+nods+1),count+1)).reshape(-1,1),1.*numpy.ones((nods,1)))))) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T)) + bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T)) if i: - bamg_geometry.SubDomains=numpy.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) + bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) + + # bamg_geometry.Vertices=np.hstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods].reshape(-1),domaini['y'][0:nods].reshape(-1),np.ones((nods)))))) + # bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.hstack((np.arange(count+1,count+nods+1).reshape(-1,),np.hstack((np.arange(count+2,count+nods+1),count+1)).reshape(-1,),1.*np.ones((nods,)))))) + # if i: + # bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) #update counter @@ -120,13 +125,13 @@ #detect whether all points of the rift are inside the domain flags=ContourToNodes(rifti['x'],rifti['y'],domain[0],0)[0] - if numpy.all(numpy.logical_not(flags)): + if np.all(np.logical_not(flags)): raise RuntimeError("one rift has all its points outside of the domain outline") - elif numpy.any(numpy.logical_not(flags)): + elif np.any(np.logical_not(flags)): #We LOTS of work to do print "Rift tip outside of or on the domain has been detected and is being processed..." #check that only one point is outside (for now) - if numpy.sum(numpy.logical_not(flags).astype(int))!=1: + if np.sum(np.logical_not(flags).astype(int))!=1: raise RuntimeError("bamg error message: only one point outside of the domain is supported yet") @@ -136,6 +141,6 @@ pass elif not flags[-1]: - rifti['x']=numpy.flipud(rifti['x']) - rifti['y']=numpy.flipud(rifti['y']) + rifti['x']=np.flipud(rifti['x']) + rifti['y']=np.flipud(rifti['y']) else: raise RuntimeError("bamg error message: only a rift tip can be outside of the domain") @@ -146,6 +151,6 @@ x2=rifti['x'][1] y2=rifti['y'][1] - for j in xrange(0,numpy.size(domain[0]['x'])-1): - if SegIntersect(numpy.array([[x1,y1],[x2,y2]]),numpy.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])): + for j in xrange(0,np.size(domain[0]['x'])-1): + if SegIntersect(np.array([[x1,y1],[x2,y2]]),np.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])): #Get position of the two nodes of the edge in domain @@ -161,11 +166,11 @@ # x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); # y=det([det([x1 y1; x2 y2]) y1-y2;det([x3 y3; x4 y4]) y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); - x=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])) - y=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])) + x=np.linalg.det(np.array([[np.linalg.det(np.array([[x1,y1],[x2,y2]])),x1-x2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),x3-x4]]))/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])) + y=np.linalg.det(np.array([[np.linalg.det(np.array([[x1,y1],[x2,y2]])),y1-y2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),y3-y4]]))/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])) segdis= sqrt((x4-x3)**2+(y4-y3)**2) - tipdis=numpy.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)]) - - if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0): + tipdis=np.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)]) + + if np.min(tipdis)/segdis < options.getfieldvalue('toltip',0): print "moving tip-domain intersection point" @@ -179,8 +184,8 @@ #OK, now we can add our own rift nods=rifti['nods']-1 - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) - bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ - numpy.array([[pos,count+1,(1+i)]]),\ - numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1)))))) + bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\ + np.array([[pos,count+1,(1+i)]]),\ + np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),(1+i)*np.ones((nods-1,1)))))) count+=nods @@ -189,22 +194,22 @@ else: #Add intersection point to Vertices - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.array([[x,y,1]]))) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.array([[x,y,1]]))) count+=1 #Decompose the crossing edge into 2 subedges - pos=numpy.nonzero(numpy.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0] + pos=np.nonzero(np.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0] if not pos: raise RuntimeError("bamg error message: a problem occurred...") - bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges[0:pos-1,:],\ - numpy.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\ - numpy.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\ + bamg_geometry.Edges=np.vstack((bamg_geometry.Edges[0:pos-1,:],\ + np.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\ + np.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\ bamg_geometry.Edges[pos+1:,:])) #OK, now we can add our own rift nods=rifti['nods']-1 - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) - bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ - numpy.array([[count,count+1,2]]),\ - numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1)))))) + bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\ + np.array([[count,count+1,2]]),\ + np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),(1+i)*np.ones((nods-1,1)))))) count+=nods @@ -213,6 +218,6 @@ else: nods=rifti['nods']-1 - bamg_geometry.Vertices=numpy.vstack(bamg_geometry.Vertices, numpy.hstack(rifti['x'][:],rifti['y'][:],numpy.ones((nods+1,1)))) - bamg_geometry.Edges =numpy.vstack(bamg_geometry.Edges, numpy.hstack(numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),i*numpy.ones((nods,1)))) + bamg_geometry.Vertices=np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:],rifti['y'][:],np.ones((nods+1,1)))) + bamg_geometry.Edges =np.vstack(bamg_geometry.Edges, np.hstack(np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),i*np.ones((nods,1)))) count=+nods+1 @@ -224,18 +229,18 @@ if all(isinstance(track,(str,unicode))): A=expread(track) - track=numpy.hstack((A.x.reshape(-1,1),A.y.reshape(-1,1))) + track=np.hstack((A.x.reshape(-1,),A.y.reshape(-1,))) else: track=float(track) #for some reason, it is of class "single" - if numpy.size(track,axis=1)==2: - track=numpy.hstack((track,3.*numpy.ones((size(track,axis=0),1)))) + if np.size(track,axis=1)==2: + track=np.hstack((track,3.*np.ones((size(track,axis=0),1)))) #only keep those inside flags=ContourToNodes(track[:,0],track[:,1],domainfile,0)[0] - track=track[numpy.nonzero(flags),:] + track=track[np.nonzero(flags),:] #Add all points to bamg_geometry - nods=numpy.size(track,axis=0) - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,track)) - bamg_geometry.Edges =numpy.vstack((bamg_geometry.Edges,numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),3.*numpy.ones((nods-1,1)))))) + nods=np.size(track,axis=0) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,track)) + bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),3.*np.ones((nods-1,1)))))) #update counter @@ -247,13 +252,13 @@ #recover RequiredVertices requiredvertices=options.getfieldvalue('RequiredVertices') #for some reason, it is of class "single" - if numpy.size(requiredvertices,axis=1)==2: - requiredvertices=numpy.hstack((requiredvertices,4.*numpy.ones((numpy.size(requiredvertices,axis=0),1)))) + if np.size(requiredvertices,axis=1)==2: + requiredvertices=np.hstack((requiredvertices,4.*np.ones((np.size(requiredvertices,axis=0),1)))) #only keep those inside flags=ContourToNodes(requiredvertices[:,0],requiredvertices[:,1],domainfile,0)[0] - requiredvertices=requiredvertices[numpy.nonzero(flags)[0],:] + requiredvertices=requiredvertices[np.nonzero(flags)[0],:] #Add all points to bamg_geometry - nods=numpy.size(requiredvertices,axis=0) - bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,requiredvertices)) + nods=np.size(requiredvertices,axis=0) + bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,requiredvertices)) #update counter @@ -275,6 +280,7 @@ bamg_mesh=bamgmesh(md.private.bamg['mesh'].__dict__) else: - bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1)))) - bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1)))) + bamg_mesh.Vertices=np.vstack((md.mesh.x,md.mesh.y,np.ones((md.mesh.numberofvertices)))).T + #bamg_mesh.Vertices=np.hstack((md.mesh.x.reshape(-1,),md.mesh.y.reshape(-1,),np.ones((md.mesh.numberofvertices,1)))) + bamg_mesh.Triangles=np.hstack((md.mesh.elements,np.ones((md.mesh.numberofelements,1)))) if isinstance(md.rifts.riftstruct,dict): @@ -286,19 +292,19 @@ bamg_options['coeff']=options.getfieldvalue('coeff',1.) bamg_options['cutoff']=options.getfieldvalue('cutoff',10.**-5) - bamg_options['err']=options.getfieldvalue('err',numpy.array([[0.01]])) + bamg_options['err']=options.getfieldvalue('err',np.array([[0.01]])) bamg_options['errg']=options.getfieldvalue('errg',0.1) - bamg_options['field']=options.getfieldvalue('field',numpy.empty((0,1))) + bamg_options['field']=options.getfieldvalue('field',np.empty((0,1))) bamg_options['gradation']=options.getfieldvalue('gradation',1.5) bamg_options['Hessiantype']=options.getfieldvalue('Hessiantype',0) bamg_options['hmin']=options.getfieldvalue('hmin',10.**-100) bamg_options['hmax']=options.getfieldvalue('hmax',10.**100) - bamg_options['hminVertices']=options.getfieldvalue('hminVertices',numpy.empty((0,1))) - bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',numpy.empty((0,1))) - bamg_options['hVertices']=options.getfieldvalue('hVertices',numpy.empty((0,1))) + bamg_options['hminVertices']=options.getfieldvalue('hminVertices',np.empty((0,1))) + bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',np.empty((0,1))) + bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1))) bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1) bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.) bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6) bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.) - bamg_options['metric']=options.getfieldvalue('metric',numpy.empty((0,1))) + bamg_options['metric']=options.getfieldvalue('metric',np.empty((0,1))) bamg_options['Metrictype']=options.getfieldvalue('Metrictype',0) bamg_options['nbjacobi']=options.getfieldvalue('nbjacobi',1) @@ -328,15 +334,15 @@ #Fill in rest of fields: - md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices=numpy.size(md.mesh.x) - md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) - md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) + md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity - md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0 + md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0 md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) #Check for orphan - if numpy.any(numpy.logical_not(numpy.in1d(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))): + if np.any(np.logical_not(np.in1d(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))): raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)") @@ -349,5 +355,5 @@ print "Checking Edge crossing..." i=0 - while (ii) + posedges=np.nonzero(geom.Edges>i) geom.Edges[posedges]=geom.Edges[posedges]-1 Index: /issm/trunk-jpl/src/m/mesh/meshconvert.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from collections import OrderedDict from BamgConvertMesh import BamgConvertMesh @@ -43,8 +43,8 @@ #Fill in rest of fields: - md.mesh.numberofelements = numpy.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices = numpy.size(md.mesh.x) - md.mesh.numberofedges = numpy.size(md.mesh.edges,axis=0) - md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.numberofelements = np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices = np.size(md.mesh.x) + md.mesh.numberofedges = np.size(md.mesh.edges,axis=0) + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices,bool) md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True Index: /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 21303) @@ -1,5 +1,5 @@ from MatlabFuncs import * from model import * -from numpy import * +from np.import * from pairoptions import * from mesh3dsurface import * Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from ElementsFromEdge import ElementsFromEdge import MatlabFuncs as m @@ -19,5 +19,5 @@ tips=rift.tips - outsidetips=tips[numpy.nonzero(flags[rift.tips-1])[0]] + outsidetips=tips[np.nonzero(flags[rift.tips-1])[0]] #we have found outsidetips, tips that touch the domain outline. go through them @@ -26,10 +26,10 @@ #find tip in the segments, take first segment (there should be 2) that holds tip, #and node_connected_to_tip is the other node on this segment: - tipindex=numpy.nonzero(rift.segments[:,0]==tip)[0] + tipindex=np.nonzero(rift.segments[:,0]==tip)[0] if tipindex: tipindex=tipindex[0] node_connected_to_tip=rift.segments[tipindex,1] else: - tipindex=numpy.nonzero(rift.segments[:,1]==tip)[0] + tipindex=np.nonzero(rift.segments[:,1]==tip)[0] tipindex=tipindex[0] node_connected_to_tip=rift.segments[tipindex,1] @@ -41,5 +41,5 @@ B=node_connected_to_tip - elements=numpy.empty(0,int) + elements=np.empty(0,int) while flags(B): #as long as B does not belong to the domain outline, keep looking. @@ -48,37 +48,37 @@ #rule out those we already detected already_detected=m.ismember(edgeelements,elements) - nextelement=edgeelements(numpy.nonzero(numpy.logical_not(already_detected))[0]) + nextelement=edgeelements(np.nonzero(np.logical_not(already_detected))[0]) #add new detected element to the list of elements we are looking for. - elements=numpy.concatenate((elements,nextelement)) + elements=np.concatenate((elements,nextelement)) #new B: - B=md.mesh.elements[nextelement-1,numpy.nonzero(numpy.logical_not(m.ismember(md.mesh.elements[nextelement-1,:],numpy.array([A,B]))))] + B=md.mesh.elements[nextelement-1,np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement-1,:],np.array([A,B]))))] #take the list of elements on one side of the rift that connect to the tip, #and duplicate the tip on them, so as to open the rift to the outside. - num=numpy.size(md.mesh.x)+1 - md.mesh.x=numpy.concatenate((md.mesh.x,md.mesh.x[tip])) - md.mesh.y=numpy.concatenate((md.mesh.y,md.mesh.y[tip])) + num=np.size(md.mesh.x)+1 + md.mesh.x=np.concatenate((md.mesh.x,md.mesh.x[tip])) + md.mesh.y=np.concatenate((md.mesh.y,md.mesh.y[tip])) md.mesh.numberofvertices=num #replace tip in elements newelements=md.mesh.elements[elements-1,:] - pos=numpy.nonzero(newelements==tip) + pos=np.nonzero(newelements==tip) newelements[pos]=num md.mesh.elements[elements-1,:]=newelements - rift.tips=numpy.concatenate((rift.tips,num)) + rift.tips=np.concatenate((rift.tips,num)) #deal with segments - tipsegments=numpy.nonzero(numpy.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0] + tipsegments=np.nonzero(np.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0] for segment_index in tipsegments: - pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0] + pos=np.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0] other_node=md.mesh.segments[segment_index,pos] if not isconnected(md.mesh.elements,other_node,tip): - pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0] + pos=np.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0] md.mesh.segments[segment_index,pos]=num #Fill in rest of fields: - md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices=numpy.size(md.mesh.x) - md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.vertexonboundary=np.zeros(np.size(md.mesh.x),bool) md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True md.rifts.numrifts=length(md.rifts.riftstruct) Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from TriMeshProcessRifts import TriMeshProcessRifts from ContourToMesh import ContourToMesh @@ -33,13 +33,13 @@ #Fill in rest of fields: numrifts=len(md.rifts.riftstruct) - md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices=numpy.size(md.mesh.x) - md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) + md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=np.size(md.mesh.x) + md.mesh.vertexonboundary=np.zeros(np.size(md.mesh.x),bool) md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True #get coordinates of rift tips for rift in md.rifts.riftstruct: - rift['tip1coordinates']=numpy.hstack((md.mesh.x[rift['tips'][0,0].astype(int)-1].reshape(-1,1),md.mesh.y[rift['tips'][0,0].astype(int)-1].reshape(-1,1))) - rift['tip2coordinates']=numpy.hstack((md.mesh.x[rift['tips'][0,1].astype(int)-1].reshape(-1,1),md.mesh.y[rift['tips'][0,1].astype(int)-1].reshape(-1,1))) + rift['tip1coordinates']=np.hstack((md.mesh.x[rift['tips'][0,0].astype(int)-1].reshape(-1,),md.mesh.y[rift['tips'][0,0].astype(int)-1].reshape(-1,))) + rift['tip2coordinates']=np.hstack((md.mesh.x[rift['tips'][0,1].astype(int)-1].reshape(-1,),md.mesh.y[rift['tips'][0,1].astype(int)-1].reshape(-1,))) #In case we have rifts that open up the domain outline, we need to open them: @@ -58,6 +58,6 @@ #get elements that are not correctly oriented in the correct direction: aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y) - pos=numpy.nonzero(aires<0)[0] - md.mesh.elements[pos,:]=numpy.hstack((md.mesh.elements[pos,1].reshape(-1,1),md.mesh.elements[pos,0].reshape(-1,1),md.mesh.elements[pos,2].reshape(-1,1))) + pos=np.nonzero(aires<0)[0] + md.mesh.elements[pos,:]=np.vstack((md.mesh.elements[pos,1],md.mesh.elements[pos,0],md.mesh.elements[pos,2])).T return md Index: /issm/trunk-jpl/src/m/mesh/roundmesh.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/roundmesh.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/roundmesh.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os from collections import OrderedDict @@ -20,12 +20,12 @@ #Get number of points on the circle - pointsonedge=numpy.floor((2.*numpy.pi*radius) / resolution) + pointsonedge=np.floor((2.*np.pi*radius) / resolution) #Calculate the cartesians coordinates of the points - x_list=numpy.ones(pointsonedge) - y_list=numpy.ones(pointsonedge) - theta=numpy.linspace(0.,2.*numpy.pi,num=pointsonedge,endpoint=False) - x_list=roundsigfig(radius*x_list*numpy.cos(theta),12) - y_list=roundsigfig(radius*y_list*numpy.sin(theta),12) + x_list=np.ones(pointsonedge) + y_list=np.ones(pointsonedge) + theta=np.linspace(0.,2.*np.pi,num=pointsonedge,endpoint=False) + x_list=roundsigfig(radius*x_list*np.cos(theta),12) + y_list=roundsigfig(radius*y_list*np.sin(theta),12) A=OrderedDict() A['x']=[x_list] @@ -39,5 +39,5 @@ #move the closest node to the center - pos=numpy.argmin(md.mesh.x**2+md.mesh.y**2) + pos=np.argmin(md.mesh.x**2+md.mesh.y**2) md.mesh.x[pos]=0. md.mesh.y[pos]=0. @@ -50,10 +50,10 @@ def roundsigfig(x,n): - digits=numpy.ceil(numpy.log10(numpy.abs(x))) + digits=np.ceil(np.log10(np.abs(x))) x=x/10.**digits - x=numpy.round(x,decimals=n) + x=np.round(x,decimals=n) x=x*10.**digits - pos=numpy.nonzero(numpy.isnan(x)) + pos=np.nonzero(np.isnan(x)) x[pos]=0. Index: /issm/trunk-jpl/src/m/mesh/squaremesh.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/squaremesh.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/squaremesh.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from NodeConnectivity import NodeConnectivity from ElementConnectivity import ElementConnectivity @@ -22,7 +22,7 @@ #initialization - index=numpy.zeros((nel,3),int) - x=numpy.zeros((nx*ny)) - y=numpy.zeros((nx*ny)) + index=np.zeros((nel,3),int) + x=np.zeros((nx*ny)) + y=np.zeros((nx*ny)) #create coordinates @@ -43,17 +43,17 @@ #Scale x and y - x=x/numpy.max(x)*Lx - y=y/numpy.max(y)*Ly + x=x/np.max(x)*Lx + y=y/np.max(y)*Ly #create segments - segments=numpy.zeros((2*(nx-1)+2*(ny-1),3),int) + segments=np.zeros((2*(nx-1)+2*(ny-1),3),int) #left edge: - segments[0:ny-1,:]=numpy.hstack((numpy.arange(2,ny+1).reshape(-1,1),numpy.arange(1,ny).reshape(-1,1),(2*numpy.arange(1,ny)-1).reshape(-1,1))) + segments[0:ny-1,:]=np.hstack((np.arange(2,ny+1).reshape(-1,),np.arange(1,ny).reshape(-1,),(2*np.arange(1,ny)-1).reshape(-1,))) #right edge: - segments[ny-1:2*(ny-1),:]=numpy.hstack((numpy.arange(ny*(nx-1)+1,nx*ny).reshape(-1,1),numpy.arange(ny*(nx-1)+2,nx*ny+1).reshape(-1,1),2*numpy.arange((ny-1)*(nx-2)+1,(nx-1)*(ny-1)+1).reshape(-1,1))) + segments[ny-1:2*(ny-1),:]=np.hstack((np.arange(ny*(nx-1)+1,nx*ny).reshape(-1,),np.arange(ny*(nx-1)+2,nx*ny+1).reshape(-1,),2*np.arange((ny-1)*(nx-2)+1,(nx-1)*(ny-1)+1).reshape(-1,))) #front edge: - segments[2*(ny-1):2*(ny-1)+(nx-1),:]=numpy.hstack((numpy.arange(2*ny,ny*nx+1,ny).reshape(-1,1),numpy.arange(ny,ny*(nx-1)+1,ny).reshape(-1,1),numpy.arange(2*(ny-1),2*(nx-1)*(ny-1)+1,2*(ny-1)).reshape(-1,1))) + segments[2*(ny-1):2*(ny-1)+(nx-1),:]=np.hstack((np.arange(2*ny,ny*nx+1,ny).reshape(-1,),np.arange(ny,ny*(nx-1)+1,ny).reshape(-1,),np.arange(2*(ny-1),2*(nx-1)*(ny-1)+1,2*(ny-1)).reshape(-1,))) #back edge - segments[2*(ny-1)+(nx-1):2*(nx-1)+2*(ny-1),:]=numpy.hstack((numpy.arange(1,(nx-2)*ny+2,ny).reshape(-1,1),numpy.arange(ny+1,ny*(nx-1)+2,ny).reshape(-1,1),numpy.arange(1,2*(nx-2)*(ny-1)+2,2*(ny-1)).reshape(-1,1))) + segments[2*(ny-1)+(nx-1):2*(nx-1)+2*(ny-1),:]=np.hstack((np.arange(1,(nx-2)*ny+2,ny).reshape(-1,),np.arange(ny+1,ny*(nx-1)+2,ny).reshape(-1,),np.arange(1,2*(nx-2)*(ny-1)+2,2*(ny-1)).reshape(-1,))) #plug coordinates and nodes @@ -62,5 +62,5 @@ md.mesh.y=y md.mesh.numberofvertices=nods - md.mesh.vertexonboundary=numpy.zeros((nods),bool) + md.mesh.vertexonboundary=np.zeros((nods),bool) md.mesh.vertexonboundary[segments[:,0:2]-1]=True Index: /issm/trunk-jpl/src/m/mesh/triangle.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/triangle.py (revision 21302) +++ /issm/trunk-jpl/src/m/mesh/triangle.py (revision 21303) @@ -1,4 +1,4 @@ import os.path -import numpy +import numpy as np from mesh2d import mesh2d from TriMesh import TriMesh @@ -56,7 +56,7 @@ #Fill in rest of fields: - md.mesh.numberofelements = numpy.size(md.mesh.elements,axis=0) - md.mesh.numberofvertices = numpy.size(md.mesh.x) - md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.numberofelements = np.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices = np.size(md.mesh.x) + md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices,bool) md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 21302) +++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 21303) @@ -49,27 +49,27 @@ def ismember(a,s): - import numpy + import numpy as np - if not isinstance(s,(tuple,list,dict,numpy.ndarray)): + if not isinstance(s,(tuple,list,dict,np.ndarray)): s=[s] - if not isinstance(a,(tuple,list,dict,numpy.ndarray)): + if not isinstance(a,(tuple,list,dict,np.ndarray)): a=[a] - if not isinstance(a,numpy.ndarray): + if not isinstance(a,np.ndarray): b=[item in s for item in a] else: - if not isinstance(s,numpy.ndarray): - b=numpy.empty_like(a) + if not isinstance(s,np.ndarray): + b=np.empty_like(a) for i,item in enumerate(a.flat): b.flat[i]=item in s else: - b=numpy.in1d(a.flat,s.flat).reshape(a.shape) + b=np.in1d(a.flat,s.flat).reshape(a.shape) return b def det(a): - import numpy + import numpy as np if a.shape==(1,): @@ -83,12 +83,12 @@ def sparse(ivec,jvec,svec,m=0,n=0,nzmax=0): - import numpy + import numpy as np if not m: - m=numpy.max(ivec) + m=np.max(ivec) if not n: - n=numpy.max(jvec) + n=np.max(jvec) - a=numpy.zeros((m,n)) + a=np.zeros((m,n)) for i,j,s in zip(ivec.reshape(-1,order='F'),jvec.reshape(-1,order='F'),svec.reshape(-1,order='F')): @@ -98,9 +98,9 @@ def heaviside(x): - import numpy + import numpy as np - y=numpy.zeros_like(x) - y[numpy.nonzero(x> 0.)]=1. - y[numpy.nonzero(x==0.)]=0.5 + y=np.zeros_like(x) + y[np.nonzero(x> 0.)]=1. + y[np.nonzero(x==0.)]=0.5 return y Index: /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 21302) +++ /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 21303) @@ -1,5 +1,6 @@ +import numpy as np + def logical_and_n(*arg): - from numpy import logical_and - + if len(arg): result=arg[0] @@ -12,6 +13,5 @@ def logical_or_n(*arg): - from numpy import logical_or - + if len(arg): result=arg[0] Index: /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py (revision 21302) +++ /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py (revision 21303) @@ -1,4 +1,4 @@ #Module import -import numpy +import numpy as np from math import isnan import MatlabFuncs as m @@ -33,5 +33,5 @@ #matrix - elif isinstance(field,numpy.ndarray): + elif isinstance(field,np.ndarray): string=displayunit(offset,name,str(field.shape),comment) Index: /issm/trunk-jpl/src/m/miscellaneous/isnans.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/isnans.py (revision 21302) +++ /issm/trunk-jpl/src/m/miscellaneous/isnans.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np def isnans(array): @@ -13,5 +13,5 @@ returnvalue=0 else: - returnvalue=numpy.isnan(array) + returnvalue=np.isnan(array) return returnvalue Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 21302) +++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 21303) @@ -1,4 +1,4 @@ import os.path -import numpy +import numpy as np import copy from NodeConnectivity import NodeConnectivity @@ -40,7 +40,7 @@ #Now, build the connectivity tables for this mesh. #Computing connectivity - if numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d: + if np.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d: md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices)[0] - if numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d: + if np.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and np.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d: md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity)[0] @@ -65,22 +65,22 @@ #get flag list of elements and nodes inside the contour nodein=ContourToMesh(elements,x,y,file,'node',1) - elemin=(numpy.sum(nodein(elements),axis=1)==numpy.size(elements,axis=1)) + elemin=(np.sum(nodein(elements),axis=1)==np.size(elements,axis=1)) #modify element connectivity - elemout=numpy.nonzero(numpy.logical_not(elemin))[0] + elemout=np.nonzero(np.logical_not(elemin))[0] elementconnectivity[elemout,:]=0 - elementconnectivity[numpy.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 + elementconnectivity[np.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 else: #get flag list of elements and nodes inside the contour - nodein=numpy.zeros(numberofvertices) - elemin=numpy.zeros(numberofelements) + nodein=np.zeros(numberofvertices) + elemin=np.zeros(numberofelements) - pos=numpy.nonzero(flags) + pos=np.nonzero(flags) elemin[pos]=1 nodein[elements[pos,:]-1]=1 #modify element connectivity - elemout=numpy.nonzero(numpy.logical_not(elemin))[0] + elemout=np.nonzero(np.logical_not(elemin))[0] elementconnectivity[elemout,:]=0 - elementconnectivity[numpy.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 + elementconnectivity[np.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 #Find element on boundary @@ -88,23 +88,23 @@ flag=copy.deepcopy(elementconnectivity) if len(args)==1: - flag[numpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]] - elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0) + flag[np.nonzero(flag)]=elemin[flag[np.nonzero(flag)]] + elementonboundary=np.logical_and(np.prod(flag,axis=1)==0,np.sum(flag,axis=1)>0) #Find segments on boundary - pos=numpy.nonzero(elementonboundary)[0] - num_segments=numpy.size(pos) - segments=numpy.zeros((num_segments*3,3),int) + pos=np.nonzero(elementonboundary)[0] + num_segments=np.size(pos) + segments=np.zeros((num_segments*3,3),int) count=0 for el1 in pos: - els2=elementconnectivity[el1,numpy.nonzero(elementconnectivity[el1,:])[0]]-1 - if numpy.size(els2)>1: - flag=numpy.intersect1d(numpy.intersect1d(elements[els2[0],:],elements[els2[1],:]),elements[el1,:]) + els2=elementconnectivity[el1,np.nonzero(elementconnectivity[el1,:])[0]]-1 + if np.size(els2)>1: + flag=np.intersect1d(np.intersect1d(elements[els2[0],:],elements[els2[1],:]),elements[el1,:]) nods1=elements[el1,:] - nods1=numpy.delete(nods1,numpy.nonzero(nods1==flag)) + nods1=np.delete(nods1,np.nonzero(nods1==flag)) segments[count,:]=[nods1[0],nods1[1],el1+1] - ord1=numpy.nonzero(nods1[0]==elements[el1,:])[0][0] - ord2=numpy.nonzero(nods1[1]==elements[el1,:])[0][0] + ord1=np.nonzero(nods1[0]==elements[el1,:])[0][0] + ord2=np.nonzero(nods1[1]==elements[el1,:])[0][0] #swap segment nodes if necessary @@ -113,20 +113,20 @@ segments[count,0]=segments[count,1] segments[count,1]=temp - segments[count,0:2]=numpy.flipud(segments[count,0:2]) + segments[count,0:2]=np.flipud(segments[count,0:2]) count+=1 else: nods1=elements[el1,:] - flag=numpy.setdiff1d(nods1,elements[els2,:]) + flag=np.setdiff1d(nods1,elements[els2,:]) for j in xrange(0,3): - nods=numpy.delete(nods1,j) - if numpy.any(m.ismember(flag,nods)): + nods=np.delete(nods1,j) + if np.any(m.ismember(flag,nods)): segments[count,:]=[nods[0],nods[1],el1+1] - ord1=numpy.nonzero(nods[0]==elements[el1,:])[0][0] - ord2=numpy.nonzero(nods[1]==elements[el1,:])[0][0] + ord1=np.nonzero(nods[0]==elements[el1,:])[0][0] + ord2=np.nonzero(nods[1]==elements[el1,:])[0][0] if ( (ord1==0 and ord2==1) or (ord1==1 and ord2==2) or (ord1==2 and ord2==0) ): temp=segments[count,0] segments[count,0]=segments[count,1] segments[count,1]=temp - segments[count,0:2]=numpy.flipud(segments[count,0:2]) + segments[count,0:2]=np.flipud(segments[count,0:2]) count+=1 segments=segments[0:count,:] Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 21302) +++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 21303) @@ -1,7 +1,5 @@ -import numpy +import numpy as np from model import model from pairoptions import pairoptions -import MatlabFuncs as m -import PythonFuncs as p from FlagElements import FlagElements @@ -33,9 +31,9 @@ #process options options=pairoptions(*args) -# options=deleteduplicates(options,1); + # options=deleteduplicates(options,1); #Find_out what kind of coupling to use coupling_method=options.getfieldvalue('coupling','tiling') - if not m.strcmpi(coupling_method,'tiling') and not m.strcmpi(coupling_method,'penalties'): + if not coupling_method in ['tiling','penalties']: raise TypeError("coupling type can only be: tiling or penalties") @@ -47,13 +45,11 @@ FSflag = FlagElements(md,options.getfieldvalue('FS','')) filltype = options.getfieldvalue('fill','none') - #Flag the elements that have not been flagged as filltype - if m.strcmpi(filltype,'SIA'): - SIAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SSAflag,HOflag)))]=True - elif m.strcmpi(filltype,'SSA'): - SSAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,HOflag,FSflag)))]=True - elif m.strcmpi(filltype,'HO'): - HOflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,SSAflag,FSflag)))]=True - + if 'SIA' in filltype: + SIAflag= ~SSAflag & ~HOflag + elif 'SSA' in filltype: + SSAflag=~SIAflag & ~HOflag & ~FSflag + elif 'HO' in filltype: + HOflag=~SIAflag & ~SSAflag & ~FSflag #check that each element has at least one flag if not any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag): @@ -63,7 +59,16 @@ if any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag>1): print "setflowequation warning message: some elements have several types, higher order type is used for them" - SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,SSAflag))]=False - SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,HOflag))]=False - SSAflag[numpy.nonzero(numpy.logical_and(SSAflag,HOflag))]=False + SIAflag[np.where(np.logical_and(SIAflag,SSAflag))]=False + SIAflag[np.where(np.logical_and(SIAflag,HOflag))]=False + SSAflag[np.where(np.logical_and(SSAflag,HOflag))]=False + + #check that L1L2 is not coupled to any other model for now + if any(L1L2flag) and any(SIAflag+SSAflag+HOflag+FSflag): + raise TypeError('L1L2 cannot be coupled to any other model') + + #Check that no HO or FS for 2d mesh + if domaintype(md.mesh)=='2Dhorizontal': + if any(FSflag+HOflag): + raise TypeError('FS and HO elements not allowed in 2d mesh, extrude it first') #FS can only be used alone for now: @@ -72,81 +77,76 @@ #Initialize node fields - nodeonSIA=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonSIA[md.mesh.elements[numpy.nonzero(SIAflag),:]-1]=True - nodeonSSA=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True - nodeonL1L2=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonL1L2[md.mesh.elements[numpy.nonzero(L1L2flag),:]-1]=True - nodeonHO=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True - nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool) - noneflag=numpy.zeros(md.mesh.numberofelements,bool) + nodeonSIA=np.zeros(md.mesh.numberofvertices,bool) + nodeonSIA[md.mesh.elements[np.where(SIAflag),:]-1]=True + nodeonSSA=np.zeros(md.mesh.numberofvertices,bool) + nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True + nodeonL1L2=np.zeros(md.mesh.numberofvertices,bool) + nodeonL1L2[md.mesh.elements[np.where(L1L2flag),:]-1]=True + nodeonHO=np.zeros(md.mesh.numberofvertices,bool) + nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True + nodeonFS=np.zeros(md.mesh.numberofvertices,bool) + noneflag=np.zeros(md.mesh.numberofelements,bool) #First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA) if any(FSflag): -# fullspcnodes=double((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy)+~isnan(md.stressbalance.spcvz))==3 | (nodeonHO & nodeonFS)); %find all the nodes on the boundary of the domain without icefront - fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)).astype(int)+ \ - numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)).astype(int)+ \ - numpy.logical_not(numpy.isnan(md.stressbalance.spcvz)).astype(int)==3, \ - numpy.logical_and(nodeonHO,nodeonFS)).astype(int) #find all the nodes on the boundary of the domain without icefront -# fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront - fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront - FSflag[numpy.nonzero(fullspcelems.reshape(-1))]=False - nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True + fullspcnodes=np.logical_or(~np.isnan(md.stressbalance.spcvx)+~np.isnan(md.stressbalance.spcvy)+~np.isnan(md.stressbalance.spcvz),np.logical_and(nodeonHO,nodeonFS)) #find all the nodes on the boundary of the domain without icefront + fullspcelems=np.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6 #find all the nodes on the boundary of the domain without icefront + FSflag[np.where(fullspcelems.reshape(-1))]=False + nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True #Then complete with NoneApproximation or the other model used if there is no FS if any(FSflag): if any(HOflag): #fill with HO - HOflag[numpy.logical_not(FSflag)]=True - nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True + HOflag[~FSflag]=True + nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True elif any(SSAflag): #fill with SSA - SSAflag[numpy.logical_not(FSflag)]=True - nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True + SSAflag[~FSflag]=True + nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True else: #fill with none - noneflag[numpy.nonzero(numpy.logical_not(FSflag))]=True + noneflag[np.where(~FSflag)]=True #Now take care of the coupling between SSA and HO - md.stressbalance.vertex_pairing=numpy.array([]) - nodeonSSAHO=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonHOFS=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonSSAFS=numpy.zeros(md.mesh.numberofvertices,bool) - SSAHOflag=numpy.zeros(md.mesh.numberofelements,bool) - SSAFSflag=numpy.zeros(md.mesh.numberofelements,bool) - HOFSflag=numpy.zeros(md.mesh.numberofelements,bool) - if m.strcmpi(coupling_method,'penalties'): + md.stressbalance.vertex_pairing=np.array([]) + nodeonSSAHO=np.zeros(md.mesh.numberofvertices,bool) + nodeonHOFS=np.zeros(md.mesh.numberofvertices,bool) + nodeonSSAFS=np.zeros(md.mesh.numberofvertices,bool) + SSAHOflag=np.zeros(md.mesh.numberofelements,bool) + SSAFSflag=np.zeros(md.mesh.numberofelements,bool) + HOFSflag=np.zeros(md.mesh.numberofelements,bool) + if coupling_method=='penalties': #Create the border nodes between HO and SSA and extrude them numnodes2d=md.mesh.numberofvertices2d numlayers=md.mesh.numberoflayers - bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1 #Nodes connected to two different types of elements + bordernodes2d=np.where(np.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1 #Nodes connected to two different types of elements #initialize and fill in penalties structure - if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))): - penalties=numpy.zeros((0,2)) + if np.all(np.logical_not(np.isnan(bordernodes2d))): + penalties=np.zeros((0,2)) for i in xrange(1,numlayers): - penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d.reshape(-1,1),bordernodes2d.reshape(-1,1)+md.mesh.numberofvertices2d*(i))))) + penalties=np.vstack((penalties,np.vstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i))).T)) md.stressbalance.vertex_pairing=penalties - elif m.strcmpi(coupling_method,'tiling'): - if any(SSAflag) and any(HOflag): #coupling SSA HO + elif coupling_method=='tiling': + if any(SSAflag) and any(HOflag): #coupling SSA HO #Find node at the border - nodeonSSAHO[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonHO))]=True + nodeonSSAHO[np.where(np.logical_and(nodeonSSA,nodeonHO))]=True #SSA elements in contact with this layer become SSAHO elements - matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAHO)[0]) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA - SSAflag[numpy.nonzero(commonelements)]=False #these elements are now SSAHOelements - SSAHOflag[numpy.nonzero(commonelements)]=True + matrixelements=nodeonSSAHO[md.mesh.elements-1] + commonelements=np.sum(matrixelements,axis=1)!=0 + commonelements[np.where(HOflag)]=False #only one layer: the elements previously in SSA + SSAflag[np.where(commonelements)]=False #these elements are now SSAHOelements + SSAHOflag[np.where(commonelements)]=True nodeonSSA[:]=False - nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True + nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(SSAHOflag)[0] - elist=numpy.zeros(numpy.size(pos),dtype=int) - elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) - elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) - pos1=numpy.nonzero(elist==1)[0] + pos=np.where(SSAHOflag)[0] + elist=np.zeros(np.size(pos),dtype=int) + elist = elist + np.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - np.sum(nodeonHO[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) + pos1=np.where(elist==1)[0] SSAflag[pos[pos1]]=True SSAHOflag[pos[pos1]]=False - pos2=numpy.nonzero(elist==-1)[0] + pos2=np.where(elist==-1)[0] HOflag[pos[pos2]]=True SSAHOflag[pos[pos2]]=False @@ -154,31 +154,31 @@ #Recompute nodes associated to these elements nodeonSSA[:]=False - nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True + nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True nodeonHO[:]=False - nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True + nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True nodeonSSAHO[:]=False - nodeonSSAHO[md.mesh.elements[numpy.nonzero(SSAHOflag),:]-1]=True + nodeonSSAHO[md.mesh.elements[np.where(SSAHOflag),:]-1]=True elif any(HOflag) and any(FSflag): #coupling HO FS #Find node at the border - nodeonHOFS[numpy.nonzero(numpy.logical_and(nodeonHO,nodeonFS))]=True + nodeonHOFS[np.where(np.logical_and(nodeonHO,nodeonFS))]=True #FS elements in contact with this layer become HOFS elements - matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonHOFS)[0]) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA - FSflag[numpy.nonzero(commonelements)]=False #these elements are now SSAHOelements - HOFSflag[numpy.nonzero(commonelements)]=True - nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True + matrixelements=nodeonHOFS[md.mesh.elements-1] + commonelements=np.sum(matrixelements,axis=1)!=0 + commonelements[np.where(HOflag)]=False #only one layer: the elements previously in SSA + FSflag[np.where(commonelements)]=False #these elements are now SSAHOelements + HOFSflag[np.where(commonelements)]=True + nodeonFS=np.zeros(md.mesh.numberofvertices,bool) + nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(HOFSflag)[0] - elist=numpy.zeros(numpy.size(pos),dtype=int) - elist = elist + numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool) - elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool) - pos1=numpy.nonzero(elist==1)[0] + pos=np.where(HOFSflag)[0] + elist=np.zeros(np.size(pos),dtype=int) + elist = elist + np.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - np.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + pos1=np.where(elist==1)[0] FSflag[pos[pos1]]=True HOFSflag[pos[pos1]]=False - pos2=numpy.nonzero(elist==-1)[0] + pos2=np.where(elist==-1)[0] HOflag[pos[pos2]]=True HOFSflag[pos[pos2]]=False @@ -186,31 +186,30 @@ #Recompute nodes associated to these elements nodeonFS[:]=False - nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True + nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True nodeonHO[:]=False - nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True + nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True nodeonHOFS[:]=False - nodeonHOFS[md.mesh.elements[numpy.nonzero(HOFSflag),:]-1]=True - + nodeonHOFS[md.mesh.elements[np.where(HOFSflag),:]-1]=True elif any(FSflag) and any(SSAflag): #Find node at the border - nodeonSSAFS[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonFS))]=True + nodeonSSAFS[np.where(np.logical_and(nodeonSSA,nodeonFS))]=True #FS elements in contact with this layer become SSAFS elements - matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAFS)[0]) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(SSAflag)]=False #only one layer: the elements previously in SSA - FSflag[numpy.nonzero(commonelements)]=False #these elements are now SSASSAelements - SSAFSflag[numpy.nonzero(commonelements)]=True - nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool) - nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True + matrixelements=nodeonSSAFS[md.mesh.elements-1] + commonelements=np.sum(matrixelements,axis=1)!=0 + commonelements[np.where(SSAflag)]=False #only one layer: the elements previously in SSA + FSflag[np.where(commonelements)]=False #these elements are now SSASSAelements + SSAFSflag[np.where(commonelements)]=True + nodeonFS=np.zeros(md.mesh.numberofvertices,bool) + nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(SSAFSflag)[0] - elist=numpy.zeros(numpy.size(pos),dtype=int) - elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) - elist = elist - numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) - pos1=numpy.nonzero(elist==1)[0] + pos=np.where(SSAFSflag)[0] + elist=np.zeros(np.size(pos),dtype=int) + elist = elist + np.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - np.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + pos1=np.where(elist==1)[0] SSAflag[pos[pos1]]=True SSAFSflag[pos[pos1]]=False - pos2=numpy.nonzero(elist==-1)[0] + pos2=np.where(elist==-1)[0] FSflag[pos[pos2]]=True SSAFSflag[pos[pos2]]=False @@ -218,9 +217,9 @@ #Recompute nodes associated to these elements nodeonSSA[:]=False - nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True + nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True nodeonFS[:]=False - nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True + nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True nodeonSSAFS[:]=False - nodeonSSAFS[md.mesh.elements[numpy.nonzero(SSAFSflag),:]-1]=True + nodeonSSAFS[md.mesh.elements[np.where(SSAFSflag),:]-1]=True elif any(FSflag) and any(SIAflag): @@ -228,14 +227,14 @@ #Create SSAHOApproximation where needed - md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int) - md.flowequation.element_equation[numpy.nonzero(noneflag)]=0 - md.flowequation.element_equation[numpy.nonzero(SIAflag)]=1 - md.flowequation.element_equation[numpy.nonzero(SSAflag)]=2 - md.flowequation.element_equation[numpy.nonzero(L1L2flag)]=3 - md.flowequation.element_equation[numpy.nonzero(HOflag)]=4 - md.flowequation.element_equation[numpy.nonzero(FSflag)]=5 - md.flowequation.element_equation[numpy.nonzero(SSAHOflag)]=6 - md.flowequation.element_equation[numpy.nonzero(SSAFSflag)]=7 - md.flowequation.element_equation[numpy.nonzero(HOFSflag)]=8 + md.flowequation.element_equation=np.zeros(md.mesh.numberofelements,int) + md.flowequation.element_equation[np.where(noneflag)]=0 + md.flowequation.element_equation[np.where(SIAflag)]=1 + md.flowequation.element_equation[np.where(SSAflag)]=2 + md.flowequation.element_equation[np.where(L1L2flag)]=3 + md.flowequation.element_equation[np.where(HOflag)]=4 + md.flowequation.element_equation[np.where(FSflag)]=5 + md.flowequation.element_equation[np.where(SSAHOflag)]=6 + md.flowequation.element_equation[np.where(SSAFSflag)]=7 + md.flowequation.element_equation[np.where(HOFSflag)]=8 #border @@ -245,25 +244,25 @@ #Create vertices_type - md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int) - pos=numpy.nonzero(nodeonSSA) + md.flowequation.vertex_equation=np.zeros(md.mesh.numberofvertices,int) + pos=np.where(nodeonSSA) md.flowequation.vertex_equation[pos]=2 - pos=numpy.nonzero(nodeonL1L2) + pos=np.where(nodeonL1L2) md.flowequation.vertex_equation[pos]=3 - pos=numpy.nonzero(nodeonHO) + pos=np.where(nodeonHO) md.flowequation.vertex_equation[pos]=4 - pos=numpy.nonzero(nodeonFS) + pos=np.where(nodeonFS) md.flowequation.vertex_equation[pos]=5 #DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority) - pos=numpy.nonzero(nodeonSIA) + pos=np.where(nodeonSIA) md.flowequation.vertex_equation[pos]=1 if any(FSflag): - pos=numpy.nonzero(numpy.logical_not(nodeonFS)) + pos=np.where(np.logical_not(nodeonFS)) if not (any(HOflag) or any(SSAflag)): md.flowequation.vertex_equation[pos]=0 - pos=numpy.nonzero(nodeonSSAHO) + pos=np.where(nodeonSSAHO) md.flowequation.vertex_equation[pos]=6 - pos=numpy.nonzero(nodeonHOFS) + pos=np.where(nodeonHOFS) md.flowequation.vertex_equation[pos]=7 - pos=numpy.nonzero(nodeonSSAFS) + pos=np.where(nodeonSSAFS) md.flowequation.vertex_equation[pos]=8 Index: /issm/trunk-jpl/src/m/parameterization/sethydrostaticmask.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/sethydrostaticmask.py (revision 21302) +++ /issm/trunk-jpl/src/m/parameterization/sethydrostaticmask.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np import os from model import model @@ -20,5 +20,5 @@ """ - if numpy.size(md.geometry.bed,axis=0)!=md.mesh.numberofvertices or numpy.size(md.geometry.base,axis=0)!=md.mesh.numberofvertices or numpy.size(md.geometry.thickness,axis=0)!=md.mesh.numberofvertices: + if np.size(md.geometry.bed,axis=0)!=md.mesh.numberofvertices or np.size(md.geometry.base,axis=0)!=md.mesh.numberofvertices or np.size(md.geometry.thickness,axis=0)!=md.mesh.numberofvertices: raise IOError("hydrostaticmask error message: fields in md.geometry do not have the right size.") @@ -27,8 +27,8 @@ #Check consistency of geometry - if any(md.geometry.base[numpy.nonzero(md.mask.groundedice_levelset>0.)]!=md.geometry.bed[numpy.nonzero(md.mask.groundedice_levelset>0.)]): + if any(md.geometry.base[np.nonzero(md.mask.groundedice_levelset>0.)]!=md.geometry.bed[np.nonzero(md.mask.groundedice_levelset>0.)]): print "WARNING: md.geometry.bed and md.geometry.base not equal on grounded ice" - if any(md.geometry.base[numpy.nonzero(md.mask.groundedice_levelset<=0.)]\n'.format(name)) fid.write('{0}=['.format(name)) Index: /issm/trunk-jpl/src/m/plot/writejsfile.py =================================================================== --- /issm/trunk-jpl/src/m/plot/writejsfile.py (revision 21302) +++ /issm/trunk-jpl/src/m/plot/writejsfile.py (revision 21303) @@ -1,3 +1,3 @@ -import numpy +import numpy as np from writejsfield import writejsfield def writejsfile(filename,model,keyname): @@ -48,5 +48,5 @@ fid.write('result["shortlabel"]="{0}";\n'.format(results[i].shortlabel)) fid.write('result["unit"]="{0}";\n'.format(results[i].unit)) - if type(results[i].data)==numpy.float64: + if type(results[i].data)==np.float64: fid.write('result["time_range"]=[{0},{1}];\n'.format(results[i].time_range[0],results[i].time_range[1])) fid.write('results["{0}"]=result;\n'.format(i)) Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 21302) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 21303) @@ -1,4 +1,3 @@ -import numpy -import math +import numpy as np import struct import pairoptions @@ -44,7 +43,7 @@ if options.exist('scale'): scale = options.getfieldvalue('scale') - if numpy.size(data) > 1 : - if numpy.size(data,0)==timeserieslength: - data=numpy.array(data) + if np.size(data) > 1 : + if np.size(data,0)==timeserieslength: + data=np.array(data) data[0:-1,:] = scale*data[0:-1,:] else: @@ -52,6 +51,6 @@ else: data = scale*data - if numpy.size(data) > 1 : - if numpy.size(data,0)==timeserieslength: + if np.size(data) > 1 : + if np.size(data,0)==timeserieslength: yts = options.getfieldvalue('yts') data[-1,:] = yts*data[-1,:] @@ -119,10 +118,10 @@ if isinstance(data,bool): - data=numpy.array([data]) + data=np.array([data]) elif isinstance(data,(list,tuple)): - data=numpy.array(data).reshape(-1,1) - if numpy.ndim(data) == 1: - if numpy.size(data): - data=data.reshape(numpy.size(data),1) + data=np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data=data.reshape(np.size(data),) else: data=data.reshape(0,0) @@ -131,9 +130,9 @@ s=data.shape #if matrix = NaN, then do not write anything - if s[0]==1 and s[1]==1 and math.isnan(data[0][0]): + if np.ndim(data)==2 and np.product(s)==1 and np.all(np.isnan(data)): s=(0,0) #first write length of record - fid.write(struct.pack('i',4+4+8*s[0]*s[1]+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type + fid.write(struct.pack('i',4+4+8*np.product(s)+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type #write data code and matrix type: @@ -142,9 +141,15 @@ #now write matrix - fid.write(struct.pack('i',s[0])) - fid.write(struct.pack('i',s[1])) - for i in xrange(s[0]): - for j in xrange(s[1]): - fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose + if np.ndim(data)==1: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',1)) + for i in xrange(s[0]): + fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose + else: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',s[1])) + for i in xrange(s[0]): + for j in xrange(s[1]): + fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose # }}} @@ -152,10 +157,10 @@ if isinstance(data,(int,long)): - data=numpy.array([data]) + data=np.array([data]) elif isinstance(data,(list,tuple)): - data=numpy.array(data).reshape(-1,1) - if numpy.ndim(data) == 1: - if numpy.size(data): - data=data.reshape(numpy.size(data),1) + data=np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data=data.reshape(np.size(data),) else: data=data.reshape(0,0) @@ -164,9 +169,9 @@ s=data.shape #if matrix = NaN, then do not write anything - if s[0]==1 and s[1]==1 and math.isnan(data[0][0]): + if np.ndim(data)==2 and np.product(s)==1 and np.all(np.isnan(data)): s=(0,0) #first write length of record - fid.write(struct.pack('i',4+4+8*s[0]*s[1]+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type + fid.write(struct.pack('i',4+4+8*np.product(s)+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type #write data code and matrix type: @@ -175,9 +180,15 @@ #now write matrix - fid.write(struct.pack('i',s[0])) - fid.write(struct.pack('i',s[1])) - for i in xrange(s[0]): - for j in xrange(s[1]): - fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose + if np.ndim(data) == 1: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',1)) + for i in xrange(s[0]): + fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose + else: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',s[1])) + for i in xrange(s[0]): + for j in xrange(s[1]): + fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose # }}} @@ -185,10 +196,10 @@ if isinstance(data,(bool,int,long,float)): - data=numpy.array([data]) + data=np.array([data]) elif isinstance(data,(list,tuple)): - data=numpy.array(data).reshape(-1,1) - if numpy.ndim(data) == 1: - if numpy.size(data): - data=data.reshape(numpy.size(data),1) + data=np.array(data).reshape(-1,) + if np.ndim(data) == 1: + if np.size(data): + data=data.reshape(np.size(data),) else: data=data.reshape(0,0) @@ -197,9 +208,9 @@ s=data.shape #if matrix = NaN, then do not write anything - if s[0]==1 and s[1]==1 and math.isnan(data[0][0]): + if np.ndim(data)==1 and np.product(s)==1 and np.all(np.isnan(data)): s=(0,0) #first write length of record - recordlength=4+4+8*s[0]*s[1]+4+4; #2 integers (32 bits) + the double matrix + code + matrix type + recordlength=4+4+8*np.product(s)+4+4; #2 integers (32 bits) + the double matrix + code + matrix type if recordlength > 2**31 : raise ValueError('field %s cannot be marshalled because it is larger than 4^31 bytes!' % enum) @@ -212,9 +223,15 @@ #now write matrix - fid.write(struct.pack('i',s[0])) - fid.write(struct.pack('i',s[1])) - for i in xrange(s[0]): - for j in xrange(s[1]): - fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose + if np.ndim(data) == 1: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',1)) + for i in xrange(s[0]): + fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose + else: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',s[1])) + for i in xrange(s[0]): + for j in xrange(s[1]): + fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose # }}} @@ -225,15 +242,15 @@ for matrix in data: if isinstance(matrix,(bool,int,long,float)): - matrix=numpy.array([matrix]) + matrix=np.array([matrix]) elif isinstance(matrix,(list,tuple)): - matrix=numpy.array(matrix).reshape(-1,1) - if numpy.ndim(matrix) == 1: - if numpy.size(matrix): - matrix=matrix.reshape(numpy.size(matrix),1) + matrix=np.array(matrix).reshape(-1,) + if np.ndim(matrix) == 1: + if np.size(matrix): + matrix=matrix.reshape(np.size(matrix),) else: matrix=matrix.reshape(0,0) s=matrix.shape - recordlength+=4*2+s[0]*s[1]*8 #row and col of matrix + matrix of doubles + recordlength+=4*2+np.product(s)*8 #row and col of matrix + matrix of doubles #write length of record @@ -249,16 +266,22 @@ for matrix in data: if isinstance(matrix,(bool,int,long,float)): - matrix=numpy.array([matrix]) + matrix=np.array([matrix]) elif isinstance(matrix,(list,tuple)): - matrix=numpy.array(matrix).reshape(-1,1) - if numpy.ndim(matrix) == 1: - matrix=matrix.reshape(numpy.size(matrix),1) + matrix=np.array(matrix).reshape(-1,) + if np.ndim(matrix) == 1: + matrix=matrix.reshape(np.size(matrix),) s=matrix.shape - fid.write(struct.pack('i',s[0])) - fid.write(struct.pack('i',s[1])) - for i in xrange(s[0]): - for j in xrange(s[1]): - fid.write(struct.pack('d',float(matrix[i][j]))) + if np.ndim(data) == 1: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',1)) + for i in xrange(s[0]): + fid.write(struct.pack('d',float(matrix[i]))) #get to the "c" convention, hence the transpose + else: + fid.write(struct.pack('i',s[0])) + fid.write(struct.pack('i',s[1])) + for i in xrange(s[0]): + for j in xrange(s[1]): + fid.write(struct.pack('d',float(matrix[i][j]))) # }}} Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 21302) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 21303) @@ -1,4 +1,4 @@ import struct -import numpy +import numpy as np from collections import OrderedDict import results as resultsclass @@ -153,5 +153,5 @@ M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] if type==1: - field=numpy.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float) + field=np.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float) elif type==2: field=struct.unpack('%ds' % M,fid.read(M))[0][:-1] @@ -159,5 +159,5 @@ N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] # field=transpose(fread(fid,[N M],'double')); - field=numpy.zeros(shape=(M,N),dtype=float) + field=np.zeros(shape=(M,N),dtype=float) for i in xrange(M): field[i,:]=struct.unpack('%dd' % N,fid.read(N*struct.calcsize('d'))) Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 21302) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 21303) @@ -7,5 +7,4 @@ from waitonlock import waitonlock from loadresultsfromcluster import loadresultsfromcluster -import MatlabFuncs as m def solve(md,solutionstring,*args): @@ -76,5 +75,5 @@ md.private.solution=solutionstring cluster=md.cluster - if m.strcmpi(options.getfieldvalue('batch','no'),'yes'): + if options.getfieldvalue('batch','no')=='yes': batch=1 else: @@ -82,5 +81,5 @@ #check model consistency - if m.strcmpi(options.getfieldvalue('checkconsistency','yes'),'yes'): + if options.getfieldvalue('checkconsistency','yes')=='yes': print "checking model consistency" ismodelselfconsistent(md) @@ -117,5 +116,5 @@ #Stop here if batch mode - if m.strcmpi(options.getfieldvalue('batch','no'),'yes'): + if options.getfieldvalue('batch','no')=='yes': print 'batch mode requested: not launching job interactively' print 'launch solution sequence on remote cluster by hand'