Changeset 21303
- Timestamp:
- 10/21/16 08:49:08 (8 years ago)
- Location:
- issm/trunk-jpl/src/m
- Files:
-
- 118 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/archive/arch.py
r21170 r21303 1 import numpy 1 import numpy as np 2 2 import math 3 3 import struct … … 150 150 def write_vector(fid,data): # {{{ 151 151 """ 152 Procedure to write a n umpyarray to an arch file152 Procedure to write a np.array to an arch file 153 153 """ 154 154 # Make sure our vector is the correct shape. 155 155 # Reshape it into a row vector if it is not correct. 156 156 if isinstance(data,(bool,int,long,float)): 157 data=n umpy.array([data])157 data=np.array([data]) 158 158 elif isinstance(data,(list,tuple)): 159 data=n umpy.array(data).reshape(-1,1)160 161 if n umpy.ndim(data) == 1:162 if n umpy.size(data):163 data=data.reshape(n umpy.size(data),1)159 data=np.array(data).reshape(-1,) 160 161 if np.ndim(data) == 1: 162 if np.size(data): 163 data=data.reshape(np.size(data),) 164 164 else: 165 165 data=data.reshape(0,0) … … 216 216 rows=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0] 217 217 cols=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0] 218 raw_data=n umpy.zeros(shape=(rows,cols),dtype=float)218 raw_data=np.zeros(shape=(rows,cols),dtype=float) 219 219 for i in xrange(rows): 220 220 raw_data[i,:]=struct.unpack('>%dd' % cols,fid.read(cols*struct.calcsize('>d'))) … … 259 259 elif format.shape[0] == 1 and format.shape[1] == 1: 260 260 code=2 261 elif isinstance(format,(list,tuple,n umpy.ndarray)):261 elif isinstance(format,(list,tuple,np.ndarray)): 262 262 code=3 263 263 else: -
issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py
r21254 r21303 1 1 import os 2 import numpy as np2 import numpy as np 3 3 def PattynSMB(md,Tf): 4 4 """ -
issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
r19527 r21303 1 1 import os 2 import numpy 2 import numpy as np 3 3 from ContourToMesh import ContourToMesh 4 4 … … 14 14 15 15 #node on Dirichlet 16 pos=n umpy.nonzero(md.mesh.vertexonboundary)17 md.stressbalance.spcvx=float('nan')*n umpy.ones(md.mesh.numberofvertices)18 md.stressbalance.spcvy=float('nan')*n umpy.ones(md.mesh.numberofvertices)19 md.stressbalance.spcvz=float('nan')*n umpy.ones(md.mesh.numberofvertices)16 pos=np.nonzero(md.mesh.vertexonboundary) 17 md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) 18 md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) 19 md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) 20 20 md.stressbalance.spcvx[pos]=0 21 21 md.stressbalance.spcvy[pos]=0 22 22 md.stressbalance.spcvz[pos]=0 23 md.stressbalance.referential=float('nan')*n umpy.ones((md.mesh.numberofvertices,6))24 md.stressbalance.loadingforce=0*n umpy.ones((md.mesh.numberofvertices,3))23 md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) 24 md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) 25 25 26 26 #Dirichlet Values 27 if isinstance(md.inversion.vx_obs,n umpy.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:27 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: 28 28 print " boundary conditions for stressbalance model: spc set as observed velocities" 29 29 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] … … 39 39 40 40 #Deal with other boundary conditions 41 if n umpy.all(numpy.isnan(md.balancethickness.thickening_rate)):42 md.balancethickness.thickening_rate=n umpy.zeros((md.mesh.numberofvertices,1))41 if np.all(np.isnan(md.balancethickness.thickening_rate)): 42 md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) 43 43 print " no balancethickness.thickening_rate specified: values set as zero" 44 md.masstransport.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))45 md.balancethickness.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))46 md.damage.spcdamage=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))44 md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 45 md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 46 md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) 47 47 48 if isinstance(md.initialization.temperature,n umpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:49 md.thermal.spctemperature=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))48 if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: 49 md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) 50 50 if hasattr(md.mesh,'vertexonsurface'): 51 pos=n umpy.nonzero(md.mesh.vertexonsurface)[0]51 pos=np.nonzero(md.mesh.vertexonsurface)[0] 52 52 md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface 53 if not isinstance(md.basalforcings.geothermalflux,n umpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:54 md.basalforcings.geothermalflux=50.*10**-3*n umpy.ones((md.mesh.numberofvertices,1)) #50 mW/m^253 if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices: 54 md.basalforcings.geothermalflux=50.*10**-3*np.ones((md.mesh.numberofvertices)) #50 mW/m^2 55 55 else: 56 56 print " no thermal boundary conditions created: no observed temperature found" -
issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py
r20910 r21303 1 1 import os 2 import numpy 2 import numpy as np 3 3 from ContourToMesh import ContourToMesh 4 4 import MatlabFuncs as m … … 27 27 raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) 28 28 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) 29 nodeonicefront=n umpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))29 nodeonicefront=np.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) 30 30 else: 31 nodeonicefront=n umpy.zeros((md.mesh.numberofvertices),bool)31 nodeonicefront=np.zeros((md.mesh.numberofvertices),bool) 32 32 33 33 # pos=find(md.mesh.vertexonboundary & ~nodeonicefront); 34 pos=n umpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]35 md.stressbalance.spcvx=float('nan')*n umpy.ones(md.mesh.numberofvertices)36 md.stressbalance.spcvy=float('nan')*n umpy.ones(md.mesh.numberofvertices)37 md.stressbalance.spcvz=float('nan')*n umpy.ones(md.mesh.numberofvertices)38 md.stressbalance.referential=float('nan')*n umpy.ones((md.mesh.numberofvertices,6))39 md.stressbalance.loadingforce=0*n umpy.ones((md.mesh.numberofvertices,3))34 pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(nodeonicefront)))[0] 35 md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) 36 md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) 37 md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) 38 md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) 39 md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) 40 40 41 41 #Icefront position 42 pos=n umpy.nonzero(nodeonicefront)[0]42 pos=np.nonzero(nodeonicefront)[0] 43 43 md.mask.ice_levelset[pos]=0 44 44 … … 53 53 values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1] 54 54 segmentsfront=1-values 55 n umpy.sum(segmentsfront,axis=1)!=numbernodesfront56 segments=n umpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]55 np.sum(segmentsfront,axis=1)!=numbernodesfront 56 segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0] 57 57 #Find all nodes for these segments and spc them 58 58 pos=md.mesh.segments[segments,0:-1]-1 59 59 else: 60 pos=n umpy.nonzero(md.mesh.vertexonboundary)[0]60 pos=np.nonzero(md.mesh.vertexonboundary)[0] 61 61 md.stressbalance.spcvx[pos]=0 62 62 md.stressbalance.spcvy[pos]=0 … … 64 64 65 65 #Dirichlet Values 66 if isinstance(md.inversion.vx_obs,n umpy.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:66 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: 67 67 #reshape to rank-2 if necessary to match spc arrays 68 if n umpy.ndim(md.inversion.vx_obs)==1:69 md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1, 1)70 if n umpy.ndim(md.inversion.vy_obs)==1:71 md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1, 1)68 if np.ndim(md.inversion.vx_obs)==1: 69 md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,) 70 if np.ndim(md.inversion.vy_obs)==1: 71 md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,) 72 72 print " boundary conditions for stressbalance model: spc set as observed velocities" 73 73 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] … … 81 81 82 82 #Deal with other boundary conditions 83 if n umpy.all(numpy.isnan(md.balancethickness.thickening_rate)):84 md.balancethickness.thickening_rate=n umpy.zeros((md.mesh.numberofvertices,1))83 if np.all(np.isnan(md.balancethickness.thickening_rate)): 84 md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) 85 85 print " no balancethickness.thickening_rate specified: values set as zero" 86 md.masstransport.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))87 md.balancethickness.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))88 md.damage.spcdamage=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))86 md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 87 md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 88 md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) 89 89 90 if isinstance(md.initialization.temperature,n umpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:91 md.thermal.spctemperature=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))90 if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: 91 md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) 92 92 if hasattr(md.mesh,'vertexonsurface'): 93 pos=n umpy.nonzero(md.mesh.vertexonsurface)[0]93 pos=np.nonzero(md.mesh.vertexonsurface)[0] 94 94 md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface 95 if not isinstance(md.basalforcings.geothermalflux,n umpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:96 md.basalforcings.geothermalflux=n umpy.zeros((md.mesh.numberofvertices,1))95 if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: 96 md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices)) 97 97 else: 98 98 print " no thermal boundary conditions created: no observed temperature found" -
issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
r20910 r21303 1 1 import os 2 import numpy 2 import numpy as np 3 3 from ContourToMesh import ContourToMesh 4 4 import MatlabFuncs as m … … 29 29 raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile) 30 30 incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) 31 vertexonicefront=n umpy.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))31 vertexonicefront=np.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1)) 32 32 else: 33 33 #Guess where the ice front is 34 vertexonfloatingice=n umpy.zeros((md.mesh.numberofvertices,1))35 pos=n umpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]34 vertexonfloatingice=np.zeros((md.mesh.numberofvertices)) 35 pos=np.nonzero(np.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0] 36 36 vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1. 37 vertexonicefront=n umpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.)37 vertexonicefront=np.logical_and(np.reshape(md.mesh.vertexonboundary,(-1,)),vertexonfloatingice>0.) 38 38 39 39 # pos=find(md.mesh.vertexonboundary & ~vertexonicefront); 40 pos=n umpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]41 if not n umpy.size(pos):40 pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(vertexonicefront)))[0] 41 if not np.size(pos): 42 42 print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually." 43 43 44 md.stressbalance.spcvx=float('nan')*n umpy.ones(md.mesh.numberofvertices)45 md.stressbalance.spcvy=float('nan')*n umpy.ones(md.mesh.numberofvertices)46 md.stressbalance.spcvz=float('nan')*n umpy.ones(md.mesh.numberofvertices)47 md.stressbalance.referential=float('nan')*n umpy.ones((md.mesh.numberofvertices,6))48 md.stressbalance.loadingforce=0*n umpy.ones((md.mesh.numberofvertices,3))44 md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices) 45 md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices) 46 md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices) 47 md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6)) 48 md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3)) 49 49 50 50 #Position of ice front 51 pos=n umpy.nonzero(vertexonicefront)[0]51 pos=np.nonzero(vertexonicefront)[0] 52 52 md.mask.ice_levelset[pos]=0 53 53 … … 62 62 values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1] 63 63 segmentsfront=1-values 64 n umpy.sum(segmentsfront,axis=1)!=numbernodesfront65 segments=n umpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]64 np.sum(segmentsfront,axis=1)!=numbernodesfront 65 segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0] 66 66 #Find all nodes for these segments and spc them 67 67 pos=md.mesh.segments[segments,0:-1]-1 68 68 else: 69 pos=n umpy.nonzero(md.mesh.vertexonboundary)[0]69 pos=np.nonzero(md.mesh.vertexonboundary)[0] 70 70 md.stressbalance.spcvx[pos]=0 71 71 md.stressbalance.spcvy[pos]=0 … … 73 73 74 74 #Dirichlet Values 75 if isinstance(md.inversion.vx_obs,n umpy.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:75 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: 76 76 print " boundary conditions for stressbalance model: spc set as observed velocities" 77 77 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos] … … 80 80 print " boundary conditions for stressbalance model: spc set as zero" 81 81 82 md.hydrology.spcwatercolumn=n umpy.zeros((md.mesh.numberofvertices,2))83 pos=n umpy.nonzero(md.mesh.vertexonboundary)[0]82 md.hydrology.spcwatercolumn=np.zeros((md.mesh.numberofvertices,2)) 83 pos=np.nonzero(md.mesh.vertexonboundary)[0] 84 84 md.hydrology.spcwatercolumn[pos,0]=1 85 85 … … 89 89 90 90 #Deal with other boundary conditions 91 if n umpy.all(numpy.isnan(md.balancethickness.thickening_rate)):92 md.balancethickness.thickening_rate=n umpy.zeros((md.mesh.numberofvertices,1))91 if np.all(np.isnan(md.balancethickness.thickening_rate)): 92 md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices)) 93 93 print " no balancethickness.thickening_rate specified: values set as zero" 94 94 95 md.masstransport.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))96 md.balancethickness.spcthickness=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))97 md.damage.spcdamage=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))95 md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 96 md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices)) 97 md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices)) 98 98 99 if isinstance(md.initialization.temperature,n umpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:100 md.thermal.spctemperature=float('nan')*n umpy.ones((md.mesh.numberofvertices,1))99 if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: 100 md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices)) 101 101 if hasattr(md.mesh,'vertexonsurface'): 102 pos=n umpy.nonzero(md.mesh.vertexonsurface)[0]102 pos=np.nonzero(md.mesh.vertexonsurface)[0] 103 103 md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface 104 if not isinstance(md.basalforcings.geothermalflux,n umpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:105 md.basalforcings.geothermalflux=n umpy.zeros((md.mesh.numberofvertices,1))106 md.basalforcings.geothermalflux[n umpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3 #50mW/m2104 if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: 105 md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices)) 106 md.basalforcings.geothermalflux[np.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3 #50mW/m2 107 107 else: 108 108 print " no thermal boundary conditions created: no observed temperature found" -
issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py
r21069 r21303 1 1 from MatlabFuncs import * 2 2 from model import * 3 from n umpyimport *3 from np.import * 4 4 5 5 def love_numbers(value,*varargin): -
issm/trunk-jpl/src/m/classes/SMBcomponents.py
r21049 r21303 38 38 def initialize(self,md): # {{{ 39 39 40 if n umpy.all(numpy.isnan(self.accumulation)):41 self.accumulation=n umpy.zeros((md.mesh.numberofvertices,1))40 if np.all(np.isnan(self.accumulation)): 41 self.accumulation=np.zeros((md.mesh.numberofvertices)) 42 42 print " no SMB.accumulation specified: values set as zero" 43 43 44 if n umpy.all(numpy.isnan(self.runoff)):45 self.runoff=n umpy.zeros((md.mesh.numberofvertices,1))44 if np.all(np.isnan(self.runoff)): 45 self.runoff=np.zeros((md.mesh.numberofvertices)) 46 46 print " no SMB.runoff specified: values set as zero" 47 47 48 if n umpy.all(numpy.isnan(self.evaporation)):49 self.evaporation=n umpy.zeros((md.mesh.numberofvertices,1))48 if np.all(np.isnan(self.evaporation)): 49 self.evaporation=np.zeros((md.mesh.numberofvertices)) 50 50 print " no SMB.evaporation specified: values set as zero" 51 51 -
issm/trunk-jpl/src/m/classes/SMBd18opdd.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 64 64 def initialize(self,md): # {{{ 65 65 66 if n umpy.all(numpy.isnan(self.s0p)):67 self.s0p=n umpy.zeros((md.mesh.numberofvertices,1))66 if np.all(np.isnan(self.s0p)): 67 self.s0p=np.zeros((md.mesh.numberofvertices)) 68 68 print " no SMBd18opdd.s0p specified: values set as zero" 69 69 70 if n umpy.all(numpy.isnan(self.s0t)):71 self.s0t=n umpy.zeros((md.mesh.numberofvertices,1))70 if np.all(np.isnan(self.s0t)): 71 self.s0t=np.zeros((md.mesh.numberofvertices)) 72 72 print " no SMBd18opdd.s0t specified: values set as zero" 73 73 … … 90 90 if 'MasstransportAnalysis' in analyses: 91 91 md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',[1]) 92 md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])93 md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])92 md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 93 md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 94 94 md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',[1]) 95 95 md = checkfield(md,'fieldname','smb.rlapslgm','>=',0,'numel',[1]) … … 98 98 md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 99 99 md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 100 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)100 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 101 101 md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',[1]) 102 102 -
issm/trunk-jpl/src/m/classes/SMBforcing.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 33 33 def initialize(self,md): # {{{ 34 34 35 if n umpy.all(numpy.isnan(self.mass_balance)):36 self.mass_balance=n umpy.zeros((md.mesh.numberofvertices,1))35 if np.all(np.isnan(self.mass_balance)): 36 self.mass_balance=np.zeros((md.mesh.numberofvertices)) 37 37 print " no SMBforcing.mass_balance specified: values set as zero" 38 38 -
issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py
r21049 r21303 40 40 def initialize(self,md): # {{{ 41 41 42 if n umpy.all(numpy.isnan(self.accumulation)):43 self.accumulation=n umpy.zeros((md.mesh.numberofvertices,1))42 if np.all(np.isnan(self.accumulation)): 43 self.accumulation=np.zeros((md.mesh.numberofvertices)) 44 44 print " no SMB.accumulation specified: values set as zero" 45 45 46 if n umpy.all(numpy.isnan(self.evaporation)):47 self.evaporation=n umpy.zeros((md.mesh.numberofvertices,1))46 if np.all(np.isnan(self.evaporation)): 47 self.evaporation=np.zeros((md.mesh.numberofvertices)) 48 48 print " no SMB.evaporation specified: values set as zero" 49 49 50 if n umpy.all(numpy.isnan(self.melt)):51 self.melt=n umpy.zeros((md.mesh.numberofvertices,1))50 if np.all(np.isnan(self.melt)): 51 self.melt=np.zeros((md.mesh.numberofvertices)) 52 52 print " no SMB.melt specified: values set as zero" 53 53 54 if n umpy.all(numpy.isnan(self.refreeze)):55 self.refreeze=n umpy.zeros((md.mesh.numberofvertices,1))54 if np.all(np.isnan(self.refreeze)): 55 self.refreeze=np.zeros((md.mesh.numberofvertices)) 56 56 print " no SMB.refreeze specified: values set as zero" 57 57 -
issm/trunk-jpl/src/m/classes/SMBpdd.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 94 94 def initialize(self,md): # {{{ 95 95 96 if n umpy.all(numpy.isnan(self.s0p)):97 self.s0p=n umpy.zeros((md.mesh.numberofvertices,1))96 if np.all(np.isnan(self.s0p)): 97 self.s0p=np.zeros((md.mesh.numberofvertices)) 98 98 print " no SMBpdd.s0p specified: values set as zero" 99 99 100 if n umpy.all(numpy.isnan(self.s0t)):101 self.s0t=n umpy.zeros((md.mesh.numberofvertices,1))100 if np.all(np.isnan(self.s0t)): 101 self.s0t=np.zeros((md.mesh.numberofvertices)) 102 102 print " no SMBpdd.s0t specified: values set as zero" 103 103 … … 119 119 if 'MasstransportAnalysis' in analyses: 120 120 md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',[1]) 121 md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])122 md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])121 md = checkfield(md,'fieldname','smb.s0p','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 122 md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 123 123 md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',[1]) 124 124 md = checkfield(md,'fieldname','smb.rlapslgm','>=',0,'numel',[1]) … … 128 128 md = checkfield(md,'fieldname','smb.precipitation','NaN',1,'Inf',1,'timeseries',1) 129 129 elif self.isdelta18o: 130 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)131 md = checkfield(md,'fieldname','smb.delta18o_surface','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)130 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 131 md = checkfield(md,'fieldname','smb.delta18o_surface','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 132 132 md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 133 133 md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 134 134 md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 135 135 md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 136 md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)137 md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)136 md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 137 md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 138 138 elif self.ismungsm: 139 139 md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) … … 141 141 md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 142 142 md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1) 143 md = checkfield(md,'fieldname','smb.Pfac','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)144 md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)145 md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,n umpy.nan],'singletimeseries',1)143 md = checkfield(md,'fieldname','smb.Pfac','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 144 md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 145 md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1) 146 146 147 147 md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1) -
issm/trunk-jpl/src/m/classes/adinversion.py
r21049 r21303 12 12 13 13 from MatlabFuncs import * 14 from n umpyimport *14 from np.import * 15 15 16 16 # ADINVERSION class definition 17 18 17 # 19 20 18 # Usage: 21 22 19 # adinversion=adinversion(); 23 24 25 20 26 21 class adinversion: … … 75 70 76 71 77 num_controls=n umpy.numel(md.inversion.control_parameters)78 num_costfunc=n umpy.size(md.inversion.cost_functions,2)72 num_controls=np.numel(md.inversion.control_parameters) 73 num_costfunc=np.size(md.inversion.cost_functions,2) 79 74 80 75 … … 95 90 96 91 if solution=='BalancethicknessSolution': 97 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices , 1],float('Nan'),1)98 md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices , 1], float('Nan'),1)92 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],float('Nan'),1) 93 md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices], float('Nan'),1) 99 94 elif solution=='BalancethicknessSoftSolution': 100 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices , 1],float('Nan'),1)95 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],float('Nan'),1) 101 96 else: 102 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices , 1],float('Nan'),1)103 if not n umpy.strcmp(domaintype(md.mesh),'2Dvertical'):104 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices , 1],float('Nan'),1)97 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],float('Nan'),1) 98 if not np.strcmp(domaintype(md.mesh),'2Dvertical'): 99 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],float('Nan'),1) 105 100 return md 106 101 … … 164 159 165 160 #process control parameters 166 num_control_parameters = n umpy.numel(self.control_parameters);161 num_control_parameters = np.numel(self.control_parameters); 167 162 WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray'); 168 163 WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer'); 169 164 170 165 #process cost functions 171 num_cost_functions=n umpy.size(self.cost_functions,2);166 num_cost_functions=np.size(self.cost_functions,2); 172 167 data=copy.deepcopy(self.cost_functions) 173 data[n umpy.nonzero(self.cost_functions==101)] =['SurfaceAbsVelMisfit'];174 data[n umpy.nonzero(self.cost_functions==102)]=['SurfaceRelVelMisfit'];175 data[n umpy.nonzero(self.cost_functions==103)]=['SurfaceLogVelMisfit'];176 data[n umpy.nonzero(self.cost_functions==104)]=['SurfaceLogVxVyMisfit'];177 data[n umpy.nonzero(self.cost_functions==105)]=['SurfaceAverageVelMisfit'];178 data[n umpy.nonzero(self.cost_functions==201)]=['ThicknessAbsMisfit'];179 data[n umpy.nonzero(self.cost_functions==501)]=['DragCoefficientAbsGradient'];180 data[n umpy.nonzero(self.cost_functions==502)]=['RheologyBbarAbsGradient'];181 data[n umpy.nonzero(self.cost_functions==503)]=['ThicknessAbsGradient'];182 data[n umpy.nonzero(self.cost_functions==504)]=['ThicknessAlongGradient'];183 data[n umpy.nonzero(self.cost_functions==505)]=['ThicknessAcrossGradient'];184 data[n umpy.nonzero(self.cost_functions==506)]=['BalancethicknessMisfit'];185 data[n umpy.nonzero(self.cost_functions==601)]=['SurfaceAbsMisfit'];186 data[n umpy.nonzero(self.cost_functions==1001)]=['Outputdefinition1'];187 data[n umpy.nonzero(self.cost_functions==1002)]=['Outputdefinition2'];188 data[n umpy.nonzero(self.cost_functions==1003)]=['Outputdefinition3'];189 data[n umpy.nonzero(self.cost_functions==1004)]=['Outputdefinition4'];190 data[n umpy.nonzero(self.cost_functions==1005)]=['Outputdefinition5'];191 data[n umpy.nonzero(self.cost_functions==1006)]=['Outputdefinition6'];192 data[n umpy.nonzero(self.cost_functions==1007)]=['Outputdefinition7'];193 data[n umpy.nonzero(self.cost_functions==1008)]=['Outputdefinition8'];194 data[n umpy.nonzero(self.cost_functions==1009)]=['Outputdefinition8'];195 data[n umpy.nonzero(self.cost_functions==1010)]=['Outputdefinition10'];168 data[np.nonzero(self.cost_functions==101)] =['SurfaceAbsVelMisfit']; 169 data[np.nonzero(self.cost_functions==102)]=['SurfaceRelVelMisfit']; 170 data[np.nonzero(self.cost_functions==103)]=['SurfaceLogVelMisfit']; 171 data[np.nonzero(self.cost_functions==104)]=['SurfaceLogVxVyMisfit']; 172 data[np.nonzero(self.cost_functions==105)]=['SurfaceAverageVelMisfit']; 173 data[np.nonzero(self.cost_functions==201)]=['ThicknessAbsMisfit']; 174 data[np.nonzero(self.cost_functions==501)]=['DragCoefficientAbsGradient']; 175 data[np.nonzero(self.cost_functions==502)]=['RheologyBbarAbsGradient']; 176 data[np.nonzero(self.cost_functions==503)]=['ThicknessAbsGradient']; 177 data[np.nonzero(self.cost_functions==504)]=['ThicknessAlongGradient']; 178 data[np.nonzero(self.cost_functions==505)]=['ThicknessAcrossGradient']; 179 data[np.nonzero(self.cost_functions==506)]=['BalancethicknessMisfit']; 180 data[np.nonzero(self.cost_functions==601)]=['SurfaceAbsMisfit']; 181 data[np.nonzero(self.cost_functions==1001)]=['Outputdefinition1']; 182 data[np.nonzero(self.cost_functions==1002)]=['Outputdefinition2']; 183 data[np.nonzero(self.cost_functions==1003)]=['Outputdefinition3']; 184 data[np.nonzero(self.cost_functions==1004)]=['Outputdefinition4']; 185 data[np.nonzero(self.cost_functions==1005)]=['Outputdefinition5']; 186 data[np.nonzero(self.cost_functions==1006)]=['Outputdefinition6']; 187 data[np.nonzero(self.cost_functions==1007)]=['Outputdefinition7']; 188 data[np.nonzero(self.cost_functions==1008)]=['Outputdefinition8']; 189 data[np.nonzero(self.cost_functions==1009)]=['Outputdefinition8']; 190 data[np.nonzero(self.cost_functions==1010)]=['Outputdefinition10']; 196 191 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray'); 197 192 WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer'); -
issm/trunk-jpl/src/m/classes/autodiff.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from dependent import dependent 3 3 from independent import independent … … 104 104 if num_dependent_objects: 105 105 names=[] 106 types=n umpy.zeros(num_dependent_objects)107 indices=n umpy.zeros(num_dependent_objects)106 types=np.zeros(num_dependent_objects) 107 indices=np.zeros(num_dependent_objects) 108 108 109 109 for i,dep in enumerate(self.dependents): … … 122 122 if num_independent_objects: 123 123 names=[None] * num_independent_objects 124 types=n umpy.zeros(num_independent_objects)124 types=np.zeros(num_independent_objects) 125 125 126 126 for i,indep in enumerate(self.independents): … … 136 136 137 137 for indep in self.independents: 138 if not n umpy.isnan(indep.fos_forward_index):138 if not np.isnan(indep.fos_forward_index): 139 139 index+=indep.fos_forward_index 140 140 break … … 153 153 154 154 for dep in self.dependents: 155 if not n umpy.isnan(dep.fos_reverse_index):155 if not np.isnan(dep.fos_reverse_index): 156 156 index+=dep.fos_reverse_index 157 157 break -
issm/trunk-jpl/src/m/classes/bamggeom.py
r14640 r21303 1 import numpy 1 import numpy as np 2 2 3 3 class bamggeom(object): … … 10 10 11 11 def __init__(self,*args): # {{{ 12 self.Vertices=n umpy.empty((0,3))13 self.Edges=n umpy.empty((0,3))14 self.TangentAtEdges=n umpy.empty((0,4))15 self.Corners=n umpy.empty((0,1))16 self.RequiredVertices=n umpy.empty((0,1))17 self.RequiredEdges=n umpy.empty((0,1))18 self.CrackedEdges=n umpy.empty((0,0))19 self.SubDomains=n umpy.empty((0,4))12 self.Vertices=np.empty((0,3)) 13 self.Edges=np.empty((0,3)) 14 self.TangentAtEdges=np.empty((0,4)) 15 self.Corners=np.empty((0,1)) 16 self.RequiredVertices=np.empty((0,1)) 17 self.RequiredEdges=np.empty((0,1)) 18 self.CrackedEdges=np.empty((0,0)) 19 self.SubDomains=np.empty((0,4)) 20 20 21 21 if not len(args): -
issm/trunk-jpl/src/m/classes/bamgmesh.py
r14640 r21303 1 import numpy 1 import numpy as np 2 2 3 3 class bamgmesh(object): … … 10 10 11 11 def __init__(self,*args): # {{{ 12 self.Vertices=n umpy.empty((0,3))13 self.Edges=n umpy.empty((0,3))14 self.Triangles=n umpy.empty((0,0))15 self.Quadrilaterals=n umpy.empty((0,0))16 self.IssmEdges=n umpy.empty((0,0))17 self.IssmSegments=n umpy.empty((0,0))18 self.VerticesOnGeomVertex=n umpy.empty((0,0))19 self.VerticesOnGeomEdge=n umpy.empty((0,0))20 self.EdgesOnGeomEdge=n umpy.empty((0,0))21 self.SubDomains=n umpy.empty((0,4))22 self.SubDomainsFromGeom=n umpy.empty((0,0))23 self.ElementConnectivity=n umpy.empty((0,0))24 self.NodalConnectivity=n umpy.empty((0,0))25 self.NodalElementConnectivity=n umpy.empty((0,0))26 self.CrackedVertices=n umpy.empty((0,0))27 self.CrackedEdges=n umpy.empty((0,0))12 self.Vertices=np.empty((0,3)) 13 self.Edges=np.empty((0,3)) 14 self.Triangles=np.empty((0,0)) 15 self.Quadrilaterals=np.empty((0,0)) 16 self.IssmEdges=np.empty((0,0)) 17 self.IssmSegments=np.empty((0,0)) 18 self.VerticesOnGeomVertex=np.empty((0,0)) 19 self.VerticesOnGeomEdge=np.empty((0,0)) 20 self.EdgesOnGeomEdge=np.empty((0,0)) 21 self.SubDomains=np.empty((0,4)) 22 self.SubDomainsFromGeom=np.empty((0,0)) 23 self.ElementConnectivity=np.empty((0,0)) 24 self.NodalConnectivity=np.empty((0,0)) 25 self.NodalElementConnectivity=np.empty((0,0)) 26 self.CrackedVertices=np.empty((0,0)) 27 self.CrackedEdges=np.empty((0,0)) 28 28 29 29 if not len(args): -
issm/trunk-jpl/src/m/classes/basalforcings.py
r21049 r21303 3 3 from checkfield import checkfield 4 4 from WriteData import WriteData 5 import numpy 5 import numpy as np 6 6 7 7 class basalforcings(object): … … 38 38 def initialize(self,md): # {{{ 39 39 40 if n umpy.all(numpy.isnan(self.groundedice_melting_rate)):41 self.groundedice_melting_rate=n umpy.zeros((md.mesh.numberofvertices,1))40 if np.all(np.isnan(self.groundedice_melting_rate)): 41 self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) 42 42 print " no basalforcings.groundedice_melting_rate specified: values set as zero" 43 43 44 if n umpy.all(numpy.isnan(self.floatingice_melting_rate)):45 self.floatingice_melting_rate=n umpy.zeros((md.mesh.numberofvertices,1))44 if np.all(np.isnan(self.floatingice_melting_rate)): 45 self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices)) 46 46 print " no basalforcings.floatingice_melting_rate specified: values set as zero" 47 47 -
issm/trunk-jpl/src/m/classes/clusters/vilje.py
r21234 r21303 107 107 fid.write('#PBS -e %s/%s/%s.errlog \n\n' % (self.executionpath,dirname,modelname)) 108 108 fid.write('export ISSM_DIR="%s/../"\n' % self.codepath) 109 fid.write('module load intelcomp/1 3.0.1\n')110 fid.write('module load mpt/2. 06\n')111 fid.write('module load petsc/3. 4.1d\n')112 fid.write('module load parmetis/4.0. 2\n')113 fid.write('module load mumps/ 4.10.0\n')109 fid.write('module load intelcomp/17.0.0\n') 110 fid.write('module load mpt/2.14\n') 111 fid.write('module load petsc/3.7.4d\n') 112 fid.write('module load parmetis/4.0.3\n') 113 fid.write('module load mumps/5.0.2\n') 114 114 fid.write('cd %s/%s/\n\n' % (self.executionpath,dirname)) 115 115 fid.write('mpiexec_mpt -np %i %s/%s %s %s/%s %s\n' % (self.np,self.codepath,executable,str(solution),self.executionpath,dirname,modelname)) -
issm/trunk-jpl/src/m/classes/dependent.py
r21049 r21303 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 from pairoptions import pairoptions 4 4 from fielddisplay import fielddisplay … … 51 51 s+="%s\n" % fielddisplay(self,'type',"type of variable ('vertex' or 'scalar')") 52 52 53 if not n umpy.isnan(self.fos_reverse_index):53 if not np.isnan(self.fos_reverse_index): 54 54 s+="%s\n" % fielddisplay(self,'fos_reverse_index',"index for fos_reverse driver of ADOLC") 55 55 if self.exp: … … 70 70 raise RuntimeError("dependent checkconsistency error: index for segments should be >=0") 71 71 72 if not n umpy.isnan(self.fos_reverse_index):72 if not np.isnan(self.fos_reverse_index): 73 73 if not strcmpi(driver,'fos_reverse'): 74 74 raise TypeError("cannot declare a dependent with a fos_reverse_index when the driver is not fos_reverse!") -
issm/trunk-jpl/src/m/classes/esa.py
r21260 r21303 2 2 from MatlabFuncs import * 3 3 from model import * 4 from n umpyimport *4 from np.import * 5 5 from checkfield import checkfield 6 6 from WriteData import WriteData -
issm/trunk-jpl/src/m/classes/flowequation.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 import copy 3 3 from project3d import project3d … … 105 105 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2]) 106 106 elif m.strcmp(md.mesh.domaintype(),'3D'): 107 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',n umpy.arange(0,8+1))108 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',n umpy.arange(0,8+1))107 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',np.arange(0,8+1)) 108 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',np.arange(0,8+1)) 109 109 else: 110 110 raise RuntimeError('mesh type not supported yet') … … 114 114 if 'StressbalanceSIAAnalysis' in analyses: 115 115 if any(self.element_equation==1): 116 if n umpy.any(numpy.logical_and(self.vertex_equation,md.mask.groundedice_levelset)):116 if np.any(np.logical_and(self.vertex_equation,md.mask.groundedice_levelset)): 117 117 print "\n !!! Warning: SIA's model is not consistent on ice shelves !!!\n" 118 118 -
issm/trunk-jpl/src/m/classes/gia.py
r21049 r21303 47 47 return md 48 48 49 md = checkfield(md,'fieldname','gia.mantle_viscosity','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1],'>',0)50 md = checkfield(md,'fieldname','gia.lithosphere_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1],'>',0)49 md = checkfield(md,'fieldname','gia.mantle_viscosity','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) 50 md = checkfield(md,'fieldname','gia.lithosphere_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) 51 51 md = checkfield(md,'fieldname','gia.cross_section_shape','numel',[1],'values',[1,2]) 52 52 -
issm/trunk-jpl/src/m/classes/groundingline.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 38 38 39 39 if not m.strcmp(self.migration,'None'): 40 if n umpy.any(numpy.isnan(md.geometry.bed)):40 if np.any(np.isnan(md.geometry.bed)): 41 41 md.checkmessage("requesting grounding line migration, but bathymetry is absent!") 42 pos=n umpy.nonzero(md.mask.groundedice_levelset>0.)[0]43 if any(n umpy.abs(md.geometry.base[pos]-md.geometry.bed[pos])>10**-10):42 pos=np.nonzero(md.mask.groundedice_levelset>0.)[0] 43 if any(np.abs(md.geometry.base[pos]-md.geometry.bed[pos])>10**-10): 44 44 md.checkmessage("base not equal to bed on grounded ice!") 45 45 if any(md.geometry.bed - md.geometry.base > 10**-9): -
issm/trunk-jpl/src/m/classes/hydrologydc.py
r21280 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from fielddisplay import fielddisplay … … 134 134 # }}} 135 135 def initialize(self,md): # {{{ 136 if n umpy.all(numpy.isnan(self.basal_moulin_input)):137 self.basal_moulin_input=n umpy.zeros((md.mesh.numberofvertices,1))136 if np.all(np.isnan(self.basal_moulin_input)): 137 self.basal_moulin_input=np.zeros((md.mesh.numberofvertices)) 138 138 print" no hydrology.basal_moulin_input specified: values set as zero" 139 139 … … 166 166 md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0.,'numel',[1]) 167 167 md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0.,'numel',[1]) 168 md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices ,1])168 md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices]) 169 169 if self.isefficientlayer==1: 170 170 md = checkfield(md,'fieldname','hydrology.spcepl_head','Inf',1,'timeseries',1) 171 md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices ,1],'values',[0,1])171 md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices],'values',[0,1]) 172 172 md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0.,'numel',[1]) 173 173 md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0.,'numel',[1]) -
issm/trunk-jpl/src/m/classes/independent.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from pairoptions import pairoptions 3 3 from fielddisplay import fielddisplay … … 17 17 self.type = '' 18 18 self.fos_forward_index = float('NaN') 19 self.fov_forward_indices = n umpy.array([])19 self.fov_forward_indices = np.array([]) 20 20 self.nods = 0 21 21 … … 34 34 s+="%s\n" % fielddisplay(self,'name',"variable name (must match corresponding String)") 35 35 s+="%s\n" % fielddisplay(self,'type',"type of variable ('vertex' or 'scalar')") 36 if not n umpy.isnan(self.fos_forward_index):36 if not np.isnan(self.fos_forward_index): 37 37 s+="%s\n" % fielddisplay(self,'fos_forward_index',"index for fos_foward driver of ADOLC") 38 if n umpy.any(numpy.logical_not(numpy.isnan(self.fov_forward_indices))):38 if np.any(np.logical_not(np.isnan(self.fov_forward_indices))): 39 39 s+="%s\n" % fielddisplay(self,'fov_forward_indices',"indices for fov_foward driver of ADOLC") 40 40 … … 46 46 # }}} 47 47 def checkconsistency(self,md,i,solution,analyses,driver): # {{{ 48 if not n umpy.isnan(self.fos_forward_index):48 if not np.isnan(self.fos_forward_index): 49 49 if not strcmpi(driver,'fos_forward'): 50 50 raise TypeError("cannot declare an independent with a fos_forward_index when the driver is not fos_forward!") -
issm/trunk-jpl/src/m/classes/initialization.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from fielddisplay import fielddisplay … … 61 61 62 62 #Lithostatic pressure by default 63 self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) 64 #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1)) 63 # self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface[:,0]-md.mesh.z) 64 #self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,)) 65 66 if np.ndim(md.geometry.surface)==2: 67 print('Reshaping md.geometry.surface for you convenience but you should fix it in you files') 68 self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface.reshape(-1,)-md.mesh.z) 69 else: 70 self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z) 71 65 72 return self 66 73 #}}} … … 70 77 def checkconsistency(self,md,solution,analyses): # {{{ 71 78 if 'StressbalanceAnalysis' in analyses: 72 if not n umpy.any(numpy.logical_or(numpy.isnan(md.initialization.vx),numpy.isnan(md.initialization.vy))):79 if not np.any(np.logical_or(np.isnan(md.initialization.vx),np.isnan(md.initialization.vy))): 73 80 md = checkfield(md,'fieldname','initialization.vx','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 74 81 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) … … 80 87 md = checkfield(md,'fieldname','initialization.vy','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 81 88 #Triangle with zero velocity 82 if n umpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\83 n umpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)):89 if np.any(np.logical_and(np.sum(np.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ 90 np.sum(np.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): 84 91 md.checkmessage("at least one triangle has all its vertices with a zero velocity") 85 92 if 'ThermalAnalysis' in analyses: … … 98 105 if 'HydrologyDCInefficientAnalysis' in analyses: 99 106 if hasattr(md.hydrology,'hydrologydc'): 100 md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])107 md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 101 108 if 'HydrologyDCEfficientAnalysis' in analyses: 102 109 if hasattr(md.hydrology,'hydrologydc'): 103 110 if md.hydrology.isefficientlayer==1: 104 md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])105 md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])111 md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 112 md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 106 113 107 114 return md … … 124 131 if md.thermal.isenthalpy: 125 132 tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure; 126 pos = n umpy.nonzero(md.initialization.temperature > tpmp)[0]133 pos = np.nonzero(md.initialization.temperature > tpmp)[0] 127 134 enthalpy = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature); 128 enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1, 1) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,1)135 enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,) 129 136 WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy'); 130 137 -
issm/trunk-jpl/src/m/classes/inversion.py
r21244 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from fielddisplay import fielddisplay … … 76 76 self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node') 77 77 self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node') 78 if not n umpy.any(numpy.isnan(self.cost_functions_coefficients)):78 if not np.any(np.isnan(self.cost_functions_coefficients)): 79 79 self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node') 80 if not n umpy.any(numpy.isnan(self.min_parameters)):80 if not np.any(np.isnan(self.min_parameters)): 81 81 self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node') 82 if not n umpy.any(numpy.isnan(self.max_parameters)):82 if not np.any(np.isnan(self.max_parameters)): 83 83 self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node') 84 84 return self … … 98 98 #maximum number of iteration in the optimization algorithm for 99 99 #each step 100 self.maxiter_per_step=20*n umpy.ones(self.nsteps)100 self.maxiter_per_step=20*np.ones(self.nsteps) 101 101 102 102 #the inversed parameter is updated as follows: … … 105 105 #inversed parameter (10^8 for B, 50 for drag) and can be decreased 106 106 #after the first iterations 107 self.gradient_scaling=50*n umpy.ones((self.nsteps,1))107 self.gradient_scaling=50*np.ones((self.nsteps,1)) 108 108 109 109 #several responses can be used: … … 113 113 #misfit(1)/misfit(0) < self.step_threshold, we go directly to 114 114 #the next step 115 self.step_threshold=.7*n umpy.ones(self.nsteps) #30 per cent decrement115 self.step_threshold=.7*np.ones(self.nsteps) #30 per cent decrement 116 116 117 117 #cost_function_threshold is a criteria to stop the control methods. … … 128 128 return md 129 129 130 num_controls=n umpy.size(md.inversion.control_parameters)131 num_costfunc=n umpy.size(md.inversion.cost_functions)130 num_controls=np.size(md.inversion.control_parameters) 131 num_costfunc=np.size(md.inversion.cost_functions) 132 132 133 133 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1]) … … 185 185 186 186 #process cost functions 187 num_cost_functions=n umpy.size(self.cost_functions)187 num_cost_functions=np.size(self.cost_functions) 188 188 data=marshallcostfunctions(self.cost_functions) 189 189 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') -
issm/trunk-jpl/src/m/classes/linearbasalforcings.py
r21049 r21303 2 2 from checkfield import checkfield 3 3 from WriteData import WriteData 4 import numpy 4 import numpy as np 5 5 6 6 class linearbasalforcings(object): … … 51 51 def initialize(self,md): # {{{ 52 52 53 if n umpy.all(numpy.isnan(self.groundedice_melting_rate)):54 self.groundedice_melting_rate=n umpy.zeros((md.mesh.numberofvertices,1))53 if np.all(np.isnan(self.groundedice_melting_rate)): 54 self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) 55 55 print " no basalforcings.groundedice_melting_rate specified: values set as zero" 56 56 … … 92 92 yts=md.constants.yts 93 93 94 floatingice_melting_rate = n umpy.zeros((md.mesh.numberofvertices,1))95 pos=n umpy.nonzero(md.geometry.base<=md.basalforcings.deepwater_elevation)94 floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) 95 pos=np.nonzero(md.geometry.base<=md.basalforcings.deepwater_elevation) 96 96 floatingice_melting_rate[pos]=md.basalforcings.deepwater_melting_rate 97 pos=n umpy.nonzero(numpy.logical_and(md.geometry.base>md.basalforcings.deepwater_elevation,md.geometry.base<md.basalforcings.upperwater_elevation))97 pos=np.nonzero(np.logical_and(md.geometry.base>md.basalforcings.deepwater_elevation,md.geometry.base<md.basalforcings.upperwater_elevation)) 98 98 floatingice_melting_rate[pos]=md.basalforcings.deepwater_melting_rate*(md.geometry.base[pos]-md.basalforcings.upperwater_elevation)/(md.basalforcings.deepwater_elevation-md.basalforcings.upperwater_elevation) 99 99 -
issm/trunk-jpl/src/m/classes/m1qn3inversion.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from fielddisplay import fielddisplay … … 98 98 self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node') 99 99 self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node') 100 if not n umpy.any(numpy.isnan(self.cost_functions_coefficients)):100 if not np.any(np.isnan(self.cost_functions_coefficients)): 101 101 self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node') 102 if not n umpy.any(numpy.isnan(self.min_parameters)):102 if not np.any(np.isnan(self.min_parameters)): 103 103 self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node') 104 if not n umpy.any(numpy.isnan(self.max_parameters)):104 if not np.any(np.isnan(self.max_parameters)): 105 105 self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node') 106 106 return self … … 137 137 return md 138 138 139 num_controls=n umpy.size(md.inversion.control_parameters)140 num_costfunc=n umpy.size(md.inversion.cost_functions)139 num_controls=np.size(md.inversion.control_parameters) 140 num_costfunc=np.size(md.inversion.cost_functions) 141 141 142 142 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1]) … … 189 189 190 190 #process cost functions 191 num_cost_functions=n umpy.size(self.cost_functions)191 num_cost_functions=np.size(self.cost_functions) 192 192 data=marshallcostfunctions(self.cost_functions) 193 193 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') -
issm/trunk-jpl/src/m/classes/mask.py
r21155 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from project3d import project3d … … 40 40 41 41 md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices]) 42 isice=n umpy.array(md.mask.ice_levelset<=0,int)43 if n umpy.sum(isice)==0:42 isice=np.array(md.mask.ice_levelset<=0,int) 43 if np.sum(isice)==0: 44 44 raise TypeError("no ice present in the domain") 45 45 -
issm/trunk-jpl/src/m/classes/maskpsl.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 import MatlabFuncs as m 3 3 from model import * … … 11 11 # Usage: 12 12 # maskpsl=maskpsl(); 13 13 14 def __init__(self,*args): # {{{ 14 15 self.groundedice_levelset = float('NaN') … … 16 17 self.ocean_levelset = float('NaN') 17 18 self.land_levelset = float('NaN') 18 19 19 20 if not len(args): 20 21 self.setdefaultparameters() … … 24 25 def __repr__(self): # {{{ 25 26 string=' masks:' 26 27 27 string="%s\n%s"%(string,fielddisplay(self,'groundedice_levelset','is ice grounded ? grounded ice if > 0, grounding line position if = 0, floating ice if < 0')) 28 28 string="%s\n%s"%(string,fielddisplay(self,'ice_levelset','presence of ice if < 0, icefront position if = 0, no ice if > 0')) 29 29 string="%s\n%s"%(string,fielddisplay(self,'ocean_levelset','is the vertex on the ocean ? yes if = 1, no if = 0')) 30 30 string="%s\n%s"%(string,fielddisplay(self,'land_levelset','is the vertex on the land ? yes if = 1, no if = 0')) 31 32 31 return string 32 33 33 # }}} 34 34 def loadobj(self): # {{{ 35 35 # This def is directly called by matlab when a model object is 36 36 # loaded. Update old properties here 37 38 37 #2014 February 5th 39 38 if numel(self.ice_levelset)>1 and all(self.ice_levelset>=0): … … 41 40 return self 42 41 # }}} 42 43 43 def setdefaultparameters(self): # {{{ 44 44 return self 45 45 46 # }}} 47 46 48 def checkconsistency(self,md,solution,analyses): # {{{ 47 48 md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices, 1]) 49 md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices, 1]) 50 md = checkfield(md,'fieldname','mask.ocean_levelset','size',[md.mesh.numberofvertices, 1]) 51 md = checkfield(md,'fieldname','mask.land_levelset','size',[md.mesh.numberofvertices, 1]) 49 md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices]) 50 md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices]) 51 md = checkfield(md,'fieldname','mask.ocean_levelset','size',[md.mesh.numberofvertices]) 52 md = checkfield(md,'fieldname','mask.land_levelset','size',[md.mesh.numberofvertices]) 52 53 isice=(md.mask.ice_levelset<=0) 53 54 if sum(isice)==0: … … 56 57 if max(md.mask.ice_levelset)<0: 57 58 print('no ice front provided') 58 59 elements=md.mesh.elements-1; elements=elements.astype(n umpy.int32, copy=False);60 icefront=n umpy.sum(md.mask.ice_levelset[elements]==0,axis=1)59 60 elements=md.mesh.elements-1; elements=elements.astype(np.int32, copy=False); 61 icefront=np.sum(md.mask.ice_levelset[elements]==0,axis=1) 61 62 if (max(icefront)==3 & m.strcmp(md.mesh.elementtype(),'Tria')) or (max(icefront==6) & m.strcmp(md.mesh.elementtype(),'Penta')): 62 63 raise RuntimeError('At least one element has all nodes on ice front, change md.mask.ice_levelset to fix it') 63 64 64 65 return md 66 65 67 # }}} 68 66 69 def extrude(self,md): # {{{ 67 70 self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node') … … 71 74 return self 72 75 # }}} 76 73 77 def mask(*args): # {{{ 74 78 if not len(args): … … 77 81 raise RuntimeError('constructor not supported') 78 82 return self 83 79 84 # }}} 85 80 86 def marshall(self,prefix,md,fid): # {{{ 81 87 WriteData(fid,prefix,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1) … … 85 91 # }}} 86 92 def savemodeljs(self,fid,modelname): # {{{ 87 88 93 writejs1Darray(fid,[modelname, '.mask.groundedice_levelset'],self.groundedice_levelset) 89 94 writejs1Darray(fid,[modelname, '.mask.ice_levelset'],self.ice_levelset) 90 95 writejs1Darray(fid,[modelname, '.mask.ocean_levelset'],self.ocean_levelset) 91 96 writejs1Darray(fid,[modelname, '.mask.land_levelset'],self.land_levelset) 97 # }}} 92 98 93 # }}} -
issm/trunk-jpl/src/m/classes/mesh2d.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 85 85 md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 86 86 md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 87 md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',n umpy.arange(1,md.mesh.numberofvertices+1))87 md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1)) 88 88 md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3]) 89 if n umpy.any(numpy.logical_not(m.ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))):89 if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): 90 90 md.checkmessage("orphan nodes have been found. Check the mesh outline") 91 91 md = checkfield(md,'fieldname','mesh.numberofelements','>',0) … … 112 112 WriteData(fid,prefix,'object',self,'class','mesh','fieldname','x','format','DoubleMat','mattype',1) 113 113 WriteData(fid,prefix,'object',self,'class','mesh','fieldname','y','format','DoubleMat','mattype',1) 114 WriteData(fid,prefix,'name','md.mesh.z','data',n umpy.zeros(self.numberofvertices),'format','DoubleMat','mattype',1);114 WriteData(fid,prefix,'name','md.mesh.z','data',np.zeros(self.numberofvertices),'format','DoubleMat','mattype',1); 115 115 WriteData(fid,prefix,'object',self,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2) 116 116 WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofelements','format','Integer') -
issm/trunk-jpl/src/m/classes/mesh3dprisms.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import * … … 108 108 md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 109 109 md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 110 md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',n umpy.arange(1,md.mesh.numberofvertices+1))110 md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1)) 111 111 md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,6]) 112 if n umpy.any(numpy.logical_not(m.ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))):112 if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))): 113 113 md.checkmessage("orphan nodes have been found. Check the mesh3dprisms outline") 114 114 md = checkfield(md,'fieldname','mesh.numberoflayers','>=',0) -
issm/trunk-jpl/src/m/classes/mesh3dsurface.py
r21049 r21303 1 1 from MatlabFuncs import * 2 2 from model import * 3 from n umpyimport *3 from np.import * 4 4 from fielddisplay import fielddisplay 5 5 from checkfield import checkfield -
issm/trunk-jpl/src/m/classes/mismipbasalforcings.py
r21271 r21303 3 3 from checkfield import checkfield 4 4 from WriteData import WriteData 5 import numpy 5 import numpy as np 6 6 7 7 class mismipbasalforcings(object): … … 21 21 self.geothermalflux = float('NaN') 22 22 23 if n umpy.all(numpy.isnan(self.groundedice_melting_rate)):24 self.groundedice_melting_rate=n umpy.zeros(md.mesh.numberofvertices)23 if np.all(np.isnan(self.groundedice_melting_rate)): 24 self.groundedice_melting_rate=np.zeros(md.mesh.numberofvertices) 25 25 print ' no basalforcings.groundedice_melting_rate specified: values set as zero' 26 26 … … 81 81 print 'WARNING: value of yts for MISMIP+ runs different from ISSM default!' 82 82 83 floatingice_melting_rate = n umpy.zeros((md.mesh.numberofvertices,1))84 floatingice_melting_rate = md.basalforcings.meltrate_factor*n umpy.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*numpy.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0)83 floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) 84 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) 85 85 86 86 WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') -
issm/trunk-jpl/src/m/classes/model.py
r21097 r21303 1 1 #module imports {{{ 2 import numpy 2 import numpy as np 3 3 import copy 4 4 import sys … … 231 231 #get elements that are inside area 232 232 flag_elem=FlagElements(md1,area) 233 if not n umpy.any(flag_elem):233 if not np.any(flag_elem): 234 234 raise RuntimeError("extracted model is empty") 235 235 236 236 #kick out all elements with 3 dirichlets 237 spc_elem=n umpy.nonzero(numpy.logical_not(flag_elem))[0]238 spc_node=n umpy.unique(md1.mesh.elements[spc_elem,:])-1239 flag=n umpy.ones(md1.mesh.numberofvertices)237 spc_elem=np.nonzero(np.logical_not(flag_elem))[0] 238 spc_node=np.unique(md1.mesh.elements[spc_elem,:])-1 239 flag=np.ones(md1.mesh.numberofvertices) 240 240 flag[spc_node]=0 241 pos=n umpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]241 pos=np.nonzero(np.logical_not(np.sum(flag[md1.mesh.elements-1],axis=1)))[0] 242 242 flag_elem[pos]=0 243 243 244 244 #extracted elements and nodes lists 245 pos_elem=n umpy.nonzero(flag_elem)[0]246 pos_node=n umpy.unique(md1.mesh.elements[pos_elem,:])-1245 pos_elem=np.nonzero(flag_elem)[0] 246 pos_node=np.unique(md1.mesh.elements[pos_elem,:])-1 247 247 248 248 #keep track of some fields 249 249 numberofvertices1=md1.mesh.numberofvertices 250 250 numberofelements1=md1.mesh.numberofelements 251 numberofvertices2=n umpy.size(pos_node)252 numberofelements2=n umpy.size(pos_elem)253 flag_node=n umpy.zeros(numberofvertices1)251 numberofvertices2=np.size(pos_node) 252 numberofelements2=np.size(pos_elem) 253 flag_node=np.zeros(numberofvertices1) 254 254 flag_node[pos_node]=1 255 255 256 256 #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) 257 Pelem=n umpy.zeros(numberofelements1,int)258 Pelem[pos_elem]=n umpy.arange(1,numberofelements2+1)259 Pnode=n umpy.zeros(numberofvertices1,int)260 Pnode[pos_node]=n umpy.arange(1,numberofvertices2+1)257 Pelem=np.zeros(numberofelements1,int) 258 Pelem[pos_elem]=np.arange(1,numberofelements2+1) 259 Pnode=np.zeros(numberofvertices1,int) 260 Pnode[pos_node]=np.arange(1,numberofvertices2+1) 261 261 262 262 #renumber the elements (some node won't exist anymore) … … 283 283 #get field 284 284 field=getattr(md1,fieldi) 285 fieldsize=n umpy.shape(field)285 fieldsize=np.shape(field) 286 286 if hasattr(field,'__dict__') and not m.ismember(fieldi,['results'])[0]: #recursive call 287 287 object_fields=vars(field) … … 289 289 #get field 290 290 field=getattr(getattr(md1,fieldi),fieldj) 291 fieldsize=n umpy.shape(field)291 fieldsize=np.shape(field) 292 292 if len(fieldsize): 293 293 #size = number of nodes * n … … 295 295 setattr(getattr(md2,fieldi),fieldj,field[pos_node]) 296 296 elif fieldsize[0]==numberofvertices1+1: 297 setattr(getattr(md2,fieldi),fieldj,n umpy.vstack((field[pos_node],field[-1,:])))297 setattr(getattr(md2,fieldi),fieldj,np.vstack((field[pos_node],field[-1,:]))) 298 298 #size = number of elements * n 299 299 elif fieldsize[0]==numberofelements1: … … 305 305 setattr(md2,fieldi,field[pos_node]) 306 306 elif fieldsize[0]==numberofvertices1+1: 307 setattr(md2,fieldi,n umpy.hstack((field[pos_node],field[-1,:])))307 setattr(md2,fieldi,np.hstack((field[pos_node],field[-1,:]))) 308 308 #size = number of elements * n 309 309 elif fieldsize[0]==numberofelements1: … … 320 320 if md1.mesh.__class__.__name__=='mesh3dprisms': 321 321 md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node] 322 pos=n umpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]323 md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos] -1]322 pos=np.where(~np.isnan(md2.mesh.uppervertex))[0] 323 md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1] 324 324 325 325 md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node] 326 pos=n umpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0]327 md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos] -1]326 pos=np.where(~np.isnan(md2.mesh.lowervertex))[0] 327 md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1] 328 328 329 329 md2.mesh.upperelements=md1.mesh.upperelements[pos_elem] 330 pos=n umpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0]331 md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos] -1]330 pos=np.where(~np.isnan(md2.mesh.upperelements))[0] 331 md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1] 332 332 333 333 md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem] 334 pos=n umpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0]335 md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos] -1]334 pos=np.where(~np.isnan(md2.mesh.lowerelements))[0] 335 md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1] 336 336 337 337 #Initial 2d mesh 338 338 if md1.mesh.__class__.__name__=='mesh3dprisms': 339 flag_elem_2d=flag_elem[n umpy.arange(0,md1.mesh.numberofelements2d)]340 pos_elem_2d=n umpy.nonzero(flag_elem_2d)[0]341 flag_node_2d=flag_node[n umpy.arange(0,md1.mesh.numberofvertices2d)]342 pos_node_2d=n umpy.nonzero(flag_node_2d)[0]343 344 md2.mesh.numberofelements2d=n umpy.size(pos_elem_2d)345 md2.mesh.numberofvertices2d=n umpy.size(pos_node_2d)339 flag_elem_2d=flag_elem[np.arange(0,md1.mesh.numberofelements2d)] 340 pos_elem_2d=np.nonzero(flag_elem_2d)[0] 341 flag_node_2d=flag_node[np.arange(0,md1.mesh.numberofvertices2d)] 342 pos_node_2d=np.nonzero(flag_node_2d)[0] 343 344 md2.mesh.numberofelements2d=np.size(pos_elem_2d) 345 md2.mesh.numberofvertices2d=np.size(pos_node_2d) 346 346 md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:] 347 347 md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1] … … 354 354 #Edges 355 355 if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'): 356 if n umpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some numpy.nans...356 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... 357 357 #renumber first two columns 358 pos=n umpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]358 pos=np.nonzero(md2.mesh.edges[:,3]!=-1)[0] 359 359 md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1] 360 360 md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1]-1] … … 362 362 md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1] 363 363 #remove edges when the 2 vertices are not in the domain. 364 md2.mesh.edges=md2.mesh.edges[n umpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]364 md2.mesh.edges=md2.mesh.edges[np.nonzero(np.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:] 365 365 #Replace all zeros by -1 in the last two columns 366 pos=n umpy.nonzero(md2.mesh.edges[:,2]==0)[0]366 pos=np.nonzero(md2.mesh.edges[:,2]==0)[0] 367 367 md2.mesh.edges[pos,2]=-1 368 pos=n umpy.nonzero(md2.mesh.edges[:,3]==0)[0]368 pos=np.nonzero(md2.mesh.edges[:,3]==0)[0] 369 369 md2.mesh.edges[pos,3]=-1 370 370 #Invert -1 on the third column with last column (Also invert first two columns!!) 371 pos=n umpy.nonzero(md2.mesh.edges[:,2]==-1)[0]371 pos=np.nonzero(md2.mesh.edges[:,2]==-1)[0] 372 372 md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3] 373 373 md2.mesh.edges[pos,3]=-1 … … 376 376 md2.mesh.edges[pos,0]=values 377 377 #Finally remove edges that do not belong to any element 378 pos=n umpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]379 md2.mesh.edges=n umpy.delete(md2.mesh.edges,pos,axis=0)378 pos=np.nonzero(np.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0] 379 md2.mesh.edges=np.delete(md2.mesh.edges,pos,axis=0) 380 380 381 381 #Penalties 382 if n umpy.any(numpy.logical_not(numpy.isnan(md2.stressbalance.vertex_pairing))):383 for i in xrange(n umpy.size(md1.stressbalance.vertex_pairing,axis=0)):382 if np.any(np.logical_not(np.isnan(md2.stressbalance.vertex_pairing))): 383 for i in xrange(np.size(md1.stressbalance.vertex_pairing,axis=0)): 384 384 md2.stressbalance.vertex_pairing[i,:]=Pnode[md1.stressbalance.vertex_pairing[i,:]] 385 md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[n umpy.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:]386 if n umpy.any(numpy.logical_not(numpy.isnan(md2.masstransport.vertex_pairing))):387 for i in xrange(n umpy.size(md1.masstransport.vertex_pairing,axis=0)):385 md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[np.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:] 386 if np.any(np.logical_not(np.isnan(md2.masstransport.vertex_pairing))): 387 for i in xrange(np.size(md1.masstransport.vertex_pairing,axis=0)): 388 388 md2.masstransport.vertex_pairing[i,:]=Pnode[md1.masstransport.vertex_pairing[i,:]] 389 md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[n umpy.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:]389 md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[np.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:] 390 390 391 391 #recreate segments … … 394 394 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)[0] 395 395 md2.mesh.segments=contourenvelope(md2) 396 md2.mesh.vertexonboundary=n umpy.zeros(numberofvertices2,bool)396 md2.mesh.vertexonboundary=np.zeros(numberofvertices2,bool) 397 397 md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True 398 398 else: … … 401 401 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)[0] 402 402 segments=contourenvelope(md2) 403 md2.mesh.vertexonboundary=n umpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)403 md2.mesh.vertexonboundary=np.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) 404 404 md2.mesh.vertexonboundary[segments[:,0:2]-1]=True 405 md2.mesh.vertexonboundary=n umpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)405 md2.mesh.vertexonboundary=np.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) 406 406 #Then do it for 3d as usual 407 407 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)[0] … … 410 410 #Boundary conditions: Dirichlets on new boundary 411 411 #Catch the elements that have not been extracted 412 orphans_elem=n umpy.nonzero(numpy.logical_not(flag_elem))[0]413 orphans_node=n umpy.unique(md1.mesh.elements[orphans_elem,:])-1412 orphans_elem=np.nonzero(np.logical_not(flag_elem))[0] 413 orphans_node=np.unique(md1.mesh.elements[orphans_elem,:])-1 414 414 #Figure out which node are on the boundary between md2 and md1 415 nodestoflag1=n umpy.intersect1d(orphans_node,pos_node)415 nodestoflag1=np.intersect1d(orphans_node,pos_node) 416 416 nodestoflag2=Pnode[nodestoflag1].astype(int)-1 417 if n umpy.size(md1.stressbalance.spcvx)>1 and numpy.size(md1.stressbalance.spcvy)>2 and numpy.size(md1.stressbalance.spcvz)>2:418 if n umpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1:417 if np.size(md1.stressbalance.spcvx)>1 and np.size(md1.stressbalance.spcvy)>2 and np.size(md1.stressbalance.spcvz)>2: 418 if np.size(md1.inversion.vx_obs)>1 and np.size(md1.inversion.vy_obs)>1: 419 419 md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] 420 420 md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2] 421 421 else: 422 md2.stressbalance.spcvx[nodestoflag2]=n umpy.nan423 md2.stressbalance.spcvy[nodestoflag2]=n umpy.nan422 md2.stressbalance.spcvx[nodestoflag2]=np.nan 423 md2.stressbalance.spcvy[nodestoflag2]=np.nan 424 424 print "\n!! extract warning: spc values should be checked !!\n\n" 425 425 #put 0 for vz 426 426 md2.stressbalance.spcvz[nodestoflag2]=0 427 if n umpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))):428 md2.thermal.spctemperature[nodestoflag2 ,0]=1427 if np.any(np.logical_not(np.isnan(md1.thermal.spctemperature))): 428 md2.thermal.spctemperature[nodestoflag2]=1 429 429 430 430 #Results fields … … 441 441 #get subfields 442 442 for solutionsubfield,subfield in fieldi.__dict__.iteritems(): 443 if n umpy.size(subfield)==numberofvertices1:443 if np.size(subfield)==numberofvertices1: 444 444 setattr(fieldr,solutionsubfield,subfield[pos_node]) 445 elif n umpy.size(subfield)==numberofelements1:445 elif np.size(subfield)==numberofelements1: 446 446 setattr(fieldr,solutionsubfield,subfield[pos_elem]) 447 447 else: … … 455 455 #get subfields 456 456 for solutionsubfield,subfield in field.__dict__.iteritems(): 457 if n umpy.size(subfield)==numberofvertices1:457 if np.size(subfield)==numberofvertices1: 458 458 setattr(fieldr,solutionsubfield,subfield[pos_node]) 459 elif n umpy.size(subfield)==numberofelements1:459 elif np.size(subfield)==numberofelements1: 460 460 setattr(fieldr,solutionsubfield,subfield[pos_elem]) 461 461 else: … … 510 510 raise TypeError("extrusionexponent must be >=0") 511 511 numlayers=args[0] 512 extrusionlist=(n umpy.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1]512 extrusionlist=(np.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1] 513 513 514 514 elif len(args)==3: #two polynomial laws … … 520 520 raise TypeError("lower and upper extrusionexponents must be >=0") 521 521 522 lowerextrusionlist=(n umpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2.523 upperextrusionlist=(n umpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2.524 extrusionlist=n umpy.unique(numpy.concatenate((lowerextrusionlist,1.-upperextrusionlist)))522 lowerextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2. 523 upperextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2. 524 extrusionlist=np.unique(np.concatenate((lowerextrusionlist,1.-upperextrusionlist))) 525 525 526 526 if numlayers<2: … … 550 550 md.mesh.extractedelements = mesh2d.extractedelements 551 551 552 x3d=n umpy.empty((0))553 y3d=n umpy.empty((0))554 z3d=n umpy.empty((0)) #the lower node is on the bed552 x3d=np.empty((0)) 553 y3d=np.empty((0)) 554 z3d=np.empty((0)) #the lower node is on the bed 555 555 thickness3d=md.geometry.thickness #thickness and bed for these nodes 556 556 bed3d=md.geometry.base … … 558 558 #Create the new layers 559 559 for i in xrange(numlayers): 560 x3d=n umpy.concatenate((x3d,md.mesh.x))561 y3d=n umpy.concatenate((y3d,md.mesh.y))560 x3d=np.concatenate((x3d,md.mesh.x)) 561 y3d=np.concatenate((y3d,md.mesh.y)) 562 562 #nodes are distributed between bed and surface accordingly to the given exponent 563 z3d=n umpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))564 number_nodes3d=n umpy.size(x3d) #number of 3d nodes for the non extruded part of the mesh563 z3d=np.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1))) 564 number_nodes3d=np.size(x3d) #number of 3d nodes for the non extruded part of the mesh 565 565 566 566 #Extrude elements 567 elements3d=n umpy.empty((0,6),int)567 elements3d=np.empty((0,6),int) 568 568 for i in xrange(numlayers-1): 569 elements3d=n umpy.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 part570 number_el3d=n umpy.size(elements3d,axis=0) #number of 3d nodes for the non extruded part of the mesh569 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 570 number_el3d=np.size(elements3d,axis=0) #number of 3d nodes for the non extruded part of the mesh 571 571 572 572 #Keep a trace of lower and upper nodes 573 lowervertex= -1*numpy.ones(number_nodes3d,int)574 uppervertex= -1*numpy.ones(number_nodes3d,int)575 lowervertex[md.mesh.numberofvertices:]=n umpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)576 uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=n umpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)573 lowervertex=np.nan*np.ones(number_nodes3d,int) 574 uppervertex=np.nan*np.ones(number_nodes3d,int) 575 lowervertex[md.mesh.numberofvertices:]=np.arange(1,(numlayers-1)*md.mesh.numberofvertices+1) 576 uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=np.arange(md.mesh.numberofvertices+1,number_nodes3d+1) 577 577 md.mesh.lowervertex=lowervertex 578 578 md.mesh.uppervertex=uppervertex 579 579 580 580 #same for lower and upper elements 581 lowerelements= -1*numpy.ones(number_el3d,int)582 upperelements= -1*numpy.ones(number_el3d,int)583 lowerelements[md.mesh.numberofelements:]=n umpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)584 upperelements[:(numlayers-2)*md.mesh.numberofelements]=n umpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)581 lowerelements=np.nan*np.ones(number_el3d,int) 582 upperelements=np.nan*np.ones(number_el3d,int) 583 lowerelements[md.mesh.numberofelements:]=np.arange(1,(numlayers-2)*md.mesh.numberofelements+1) 584 upperelements[:(numlayers-2)*md.mesh.numberofelements]=np.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1) 585 585 md.mesh.lowerelements=lowerelements 586 586 md.mesh.upperelements=upperelements … … 605 605 606 606 #bedinfo and surface info 607 md.mesh.vertexonbase=project3d(md,'vector',n umpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)608 md.mesh.vertexonsurface=project3d(md,'vector',n umpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)607 md.mesh.vertexonbase=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1) 608 md.mesh.vertexonsurface=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers) 609 609 md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node') 610 610 … … 630 630 631 631 #connectivity 632 md.mesh.elementconnectivity=n umpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))633 md.mesh.elementconnectivity[n umpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1634 if not n umpy.isnan(md.mesh.elementconnectivity).all():632 md.mesh.elementconnectivity=np.tile(md.mesh.elementconnectivity,(numlayers-1,1)) 633 md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1 634 if not np.isnan(md.mesh.elementconnectivity).all(): 635 635 for i in xrange(1,numlayers-1): 636 636 md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \ 637 637 =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d 638 md.mesh.elementconnectivity[n umpy.nonzero(md.mesh.elementconnectivity<0)]=0638 md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity<0)]=0 639 639 640 640 md.materials.extrude(md) … … 674 674 675 675 #observations 676 if not n umpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)677 if not n umpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)678 if not n umpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)679 if not n umpy.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers)680 if isinstance(md.inversion.min_parameters,n umpy.ndarray):676 if not np.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) 677 if not np.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) 678 if not np.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) 679 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) 680 if isinstance(md.inversion.min_parameters,np.ndarray): 681 681 if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) 682 if isinstance(md.inversion.max_parameters,n umpy.ndarray):682 if isinstance(md.inversion.max_parameters,np.ndarray): 683 683 if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers) 684 if not n umpy.isnan(md.smb.mass_balance).all():684 if not np.isnan(md.smb.mass_balance).all(): 685 685 md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) 686 686 687 687 #results 688 if not n umpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)689 if not n umpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy)690 if not n umpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz)691 if not n umpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)692 if not n umpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)693 if not n umpy.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1)694 if not n umpy.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1)695 if not n umpy.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1)696 if not n umpy.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1)688 if not np.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) 689 if not np.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) 690 if not np.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) 691 if not np.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) 692 if not np.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) 693 if not np.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1) 694 if not np.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) 695 if not np.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1) 696 if not np.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1) 697 697 698 698 #gia 699 if not n umpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)700 if not n umpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)699 if not np.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) 700 if not np.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) 701 701 702 702 #elementstype 703 if not n umpy.isnan(md.flowequation.element_equation).all():703 if not np.isnan(md.flowequation.element_equation).all(): 704 704 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1) 705 705 md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1) … … 725 725 md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers) 726 726 md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers) 727 if not n umpy.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)727 if not np.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) 728 728 md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1) 729 729 … … 752 752 md.geometry.thickness=project2d(md,md.geometry.thickness,1) 753 753 md.geometry.base=project2d(md,md.geometry.base,1) 754 if isinstance(md.geometry.bed,n umpy.ndarray):754 if isinstance(md.geometry.bed,np.ndarray): 755 755 md.geometry.bed=project2d(md,md.geometry.bed,1) 756 756 md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1) … … 758 758 759 759 #lat long 760 if isinstance(md.mesh.lat,n umpy.ndarray):760 if isinstance(md.mesh.lat,np.ndarray): 761 761 if md.mesh.lat.size==md.mesh.numberofvertices: md.mesh.lat=project2d(md,md.mesh.lat,1) 762 if isinstance(md.mesh.long,n umpy.ndarray):762 if isinstance(md.mesh.long,np.ndarray): 763 763 if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1) 764 764 … … 770 770 mesh.numberofelements=md.mesh.numberofelements2d 771 771 mesh.elements=md.mesh.elements2d 772 if not n umpy.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)773 if not n umpy.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)772 if not np.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) 773 if not np.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) 774 774 md.mesh=mesh 775 775 -
issm/trunk-jpl/src/m/classes/outputdefinition.py
r21254 r21303 2 2 from checkfield import checkfield 3 3 from WriteData import WriteData 4 import numpy as np4 import numpy as np 5 5 6 6 class outputdefinition(object): -
issm/trunk-jpl/src/m/classes/qmu.py
r21254 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from collections import OrderedDict … … 105 105 # }}} 106 106 def extrude(self,md): # {{{ 107 self.partition=project3d(md,'vector',n umpy.transpose(self.partition),'type','node')107 self.partition=project3d(md,'vector',np.transpose(self.partition),'type','node') 108 108 return self 109 109 #}}} … … 120 120 md.checkmessage("concurrency should be set to 1 when running dakota in library mode") 121 121 if md.qmu.partition: 122 if not n umpy.size(md.qmu.partition)==md.mesh.numberofvertices:122 if not np.size(md.qmu.partition)==md.mesh.numberofvertices: 123 123 md.checkmessage("user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1") 124 124 if not min(md.qmu.partition)==0: -
issm/trunk-jpl/src/m/classes/results.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from pairoptions import pairoptions 3 3 from fielddisplay import fielddisplay -
issm/trunk-jpl/src/m/classes/rifts.py
r21245 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 33 33 #}}} 34 34 def checkconsistency(self,md,solution,analyses): # {{{ 35 if (not self.riftstruct) or n umpy.any(isnans(self.riftstruct)):35 if (not self.riftstruct) or np.any(isnans(self.riftstruct)): 36 36 numrifts=0 37 37 else: … … 43 43 if not isinstance(self.riftstruct,list): 44 44 md.checkmessage("rifts.riftstruct should be a structure!") 45 if n umpy.any(md.mesh.segmentmarkers>=2):45 if np.any(md.mesh.segmentmarkers>=2): 46 46 #We have segments with rift markers, but no rift structure! 47 47 md.checkmessage("model should be processed for rifts (run meshprocessrifts)!") … … 49 49 md = checkfield(md,'fieldname',"rifts.riftstruct[%d]['fill']" % i,'values',['Water','Air','Ice','Melange',0,1,2,3]) 50 50 else: 51 if self.riftstruct and n umpy.any(numpy.logical_not(isnans(self.riftstruct))):51 if self.riftstruct and np.any(np.logical_not(isnans(self.riftstruct))): 52 52 md.checkmessage("riftstruct should be NaN since numrifts is 0!") 53 53 … … 57 57 58 58 #Process rift info 59 if (not self.riftstruct) or n umpy.any(isnans(self.riftstruct)):59 if (not self.riftstruct) or np.any(isnans(self.riftstruct)): 60 60 numrifts=0 61 61 else: … … 64 64 numpairs=0 65 65 for rift in self.riftstruct: 66 numpairs+=n umpy.size(rift['penaltypairs'],axis=0)66 numpairs+=np.size(rift['penaltypairs'],axis=0) 67 67 68 68 # Convert strings in riftstruct to hard coded numbers … … 76 76 77 77 # 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. 78 data=n umpy.zeros((numpairs,12))78 data=np.zeros((numpairs,12)) 79 79 count=0 80 80 for rift in self.riftstruct: 81 numpairsforthisrift=n umpy.size(rift['penaltypairs'],0)81 numpairsforthisrift=np.size(rift['penaltypairs'],0) 82 82 data[count:count+numpairsforthisrift,0:7]=rift['penaltypairs'] 83 83 data[count:count+numpairsforthisrift,7]=rift['fill'] -
issm/trunk-jpl/src/m/classes/slr.py
r21295 r21303 2 2 from MatlabFuncs import * 3 3 from model import * 4 from numpy import * 4 import numpy as np 5 5 from checkfield import checkfield 6 6 from WriteData import WriteData … … 15 15 16 16 def __init__(self): # {{{ 17 self.deltathickness = NaN18 self.sealevel = NaN17 self.deltathickness = float('NaN') 18 self.sealevel = float('NaN') 19 19 self.maxiter = 0 20 20 self.reltol = 0 … … 60 60 61 61 #Convergence criterion: absolute, relative and residual 62 self.reltol= NaN#default62 self.reltol=float('NaN') #default 63 63 self.abstol=0.001 #1 mm of sea level rise 64 64 … … 95 95 return md 96 96 97 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements ,1])98 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices ,1])97 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) 98 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) 99 99 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1) 100 100 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1) -
issm/trunk-jpl/src/m/classes/steadystate.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from fielddisplay import fielddisplay 3 3 from checkfield import checkfield … … 54 54 md.checkmessage("for a steadystate computation, timestepping.time_step must be zero.") 55 55 56 if n umpy.isnan(md.stressbalance.reltol):56 if np.isnan(md.stressbalance.reltol): 57 57 md.checkmessage("for a steadystate computation, stressbalance.reltol (relative convergence criterion) must be defined!") 58 58 -
issm/trunk-jpl/src/m/classes/stressbalance.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 import sys 3 3 import copy … … 149 149 #singular solution 150 150 # if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2), 151 if not n umpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)),numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)))):151 if not np.any(np.logical_and(np.logical_not(np.isnan(md.stressbalance.spcvx)),np.logical_not(np.isnan(md.stressbalance.spcvy)))): 152 152 print "\n !!! Warning: no spc applied, model might not be well posed if no basal friction is applied, check for solution crash\n" 153 153 #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES 154 154 # if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6), 155 if n umpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=6)):155 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)): 156 156 md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values") 157 157 #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL 158 158 # if any(sum(isnan(md.stressbalance.referential),2)==0), 159 if n umpy.any(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)==0):160 pos=[i for i,item in enumerate(n umpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)) if item==0]161 # n umpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply162 # if n umpy.any(numpy.abs(numpy.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):159 if np.any(np.sum(np.isnan(md.stressbalance.referential),axis=1)==0): 160 pos=[i for i,item in enumerate(np.sum(np.isnan(md.stressbalance.referential),axis=1)) if item==0] 161 # np.inner (and np.dot) calculate all the dot product permutations, resulting in a full matrix multiply 162 # if np.any(np.abs(np.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon): 163 163 # md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal") 164 164 for item in md.stressbalance.referential[pos,:]: 165 if n umpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:165 if np.abs(np.inner(item[0:2],item[3:5]))>sys.float_info.epsilon: 166 166 md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal") 167 167 #CHECK THAT NO rotation specified for FS Grounded ice at base 168 168 if m.strcmp(md.mesh.domaintype(),'3D') and md.flowequation.isFS: 169 pos=n umpy.nonzero(numpy.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase))170 if n umpy.any(numpy.logical_not(numpy.isnan(md.stressbalance.referential[pos,:]))):169 pos=np.nonzero(np.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase)) 170 if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos,:]))): 171 171 md.checkmessage("no referential should be specified for basal vertices of grounded ice") 172 172 … … 195 195 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1) 196 196 197 if isinstance(self.loadingforce, (list, tuple, n umpy.ndarray)):197 if isinstance(self.loadingforce, (list, tuple, np.ndarray)): 198 198 lx=self.loadingforce[:,0]; 199 199 ly=self.loadingforce[:,1]; -
issm/trunk-jpl/src/m/classes/taoinversion.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from WriteData import WriteData … … 125 125 126 126 127 num_controls= n umpy.numel(md.inversion.control_parameters)128 num_costfunc= n umpy.size(md.inversion.cost_functions,2)127 num_controls= np.numel(md.inversion.control_parameters) 128 num_costfunc= np.size(md.inversion.cost_functions,2) 129 129 130 130 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1]) … … 155 155 156 156 if solution=='BalancethicknessSolution': 157 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices , 1],'NaN',1,'Inf',1)157 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) 158 158 elif solution=='BalancethicknessSoftSolution': 159 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices , 1],'NaN',1,'Inf',1)159 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) 160 160 else: 161 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices , 1],'NaN',1,'Inf',1)162 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices , 1],'NaN',1,'Inf',1)161 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) 162 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) 163 163 164 164 def marshall(self, md, fid): … … 188 188 189 189 #process control parameters 190 num_control_parameters = n umpy.numel(self.control_parameters)190 num_control_parameters = np.numel(self.control_parameters) 191 191 WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray') 192 192 WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer') 193 193 194 194 #process cost functions 195 num_cost_functions = n umpy.size(self.cost_functions,2)195 num_cost_functions = np.size(self.cost_functions,2) 196 196 data= marshallcostfunctions(self.cost_functions) 197 197 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') -
issm/trunk-jpl/src/m/classes/thermal.py
r21049 r21303 1 import numpy 1 import numpy as np 2 2 from project3d import project3d 3 3 from fielddisplay import fielddisplay 4 4 from checkfield import checkfield 5 5 from WriteData import WriteData 6 import MatlabFuncs as m7 6 8 7 class thermal(object): … … 44 43 #}}} 45 44 def extrude(self,md): # {{{ 46 self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',n umpy.nan)47 if isinstance(md.initialization.temperature,n umpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:48 self.spctemperature=float(' nan')*numpy.ones((md.mesh.numberofvertices,1))49 pos=n umpy.nonzero(md.mesh.vertexonsurface)[0]45 self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',np.nan) 46 if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices: 47 self.spctemperature=float('NaN')*np.ones((md.mesh.numberofvertices)) 48 pos=np.where(md.mesh.vertexonsurface)[0] 50 49 self.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface 51 50 return self … … 95 94 md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2]) 96 95 md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1) 96 md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) 97 97 98 if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3: 98 pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))) 99 replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) 100 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") 99 pos=np.where(~np.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])) 100 try: 101 spccol=np.size(md.thermal.spctemperature,1) 102 except IndexError: 103 spccol=1 104 replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)) 105 control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate 106 md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") 101 107 md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1]) 102 108 md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); 103 109 if(md.thermal.isenthalpy): 104 if n umpy.isnan(md.stressbalance.reltol):110 if np.isnan(md.stressbalance.reltol): 105 111 md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!") 106 112 md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero"); 107 md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) 113 108 114 109 115 return md -
issm/trunk-jpl/src/m/consistency/checkfield.py
r20913 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 from pairoptions import pairoptions … … 44 44 45 45 if isinstance(field,(bool,int,long,float)): 46 field=n umpy.array([field])46 field=np.array([field]) 47 47 48 48 #check empty … … 55 55 if options.exist('size'): 56 56 fieldsize=options.getfieldvalue('size') 57 if 58 if numpy.isnan(fieldsize[0]):57 if len(fieldsize) == 1: 58 if np.isnan(fieldsize[0]): 59 59 pass 60 elif not numpy.size(field,0)==fieldsize[0]: 61 md = md.checkmessage(options.getfieldvalue('message',\ 62 "field '%s' size should be %d" % (fieldname,fieldsize[0]))) 60 elif np.ndim(field)==1: 61 if not np.size(field)==fieldsize[0]: 62 md = md.checkmessage(options.getfieldvalue('message',"field {} size should be {}".format(fieldname,fieldsize[0]))) 63 else: 64 try: 65 exec("md.{}=field[:,0]".format(fieldname)) 66 print('{} had a bad dimension, we fixed it but you should check it'.format(fieldname)) 67 except IndexError: 68 md = md.checkmessage(options.getfieldvalue('message',"field {} should have {} dimension".format(fieldname,len(fieldsize)))) 63 69 elif len(fieldsize) == 2: 64 if n umpy.isnan(fieldsize[0]):65 if not n umpy.size(field,1)==fieldsize[1]:70 if np.isnan(fieldsize[0]): 71 if not np.size(field,1)==fieldsize[1]: 66 72 md = md.checkmessage(options.getfieldvalue('message',\ 67 73 "field '%s' should have %d columns" % (fieldname,fieldsize[1]))) 68 elif n umpy.isnan(fieldsize[1]):69 if not n umpy.size(field,0)==fieldsize[0]:74 elif np.isnan(fieldsize[1]): 75 if not np.size(field,0)==fieldsize[0]: 70 76 md = md.checkmessage(options.getfieldvalue('message',\ 71 77 "field '%s' should have %d lines" % (fieldname,fieldsize[0]))) 72 78 else: 73 if (not n umpy.size(field,0)==fieldsize[0]) or (not numpy.size(field,1)==fieldsize[1]):79 if (not np.size(field,0)==fieldsize[0]) or (not np.size(field,1)==fieldsize[1]): 74 80 md = md.checkmessage(options.getfieldvalue('message',\ 75 81 "field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1]))) … … 78 84 if options.exist('numel'): 79 85 fieldnumel=options.getfieldvalue('numel') 80 if n umpy.size(field) not in fieldnumel:86 if np.size(field) not in fieldnumel: 81 87 if len(fieldnumel)==1: 82 88 md = md.checkmessage(options.getfieldvalue('message',\ … … 91 97 #check NaN 92 98 if options.getfieldvalue('NaN',0): 93 if n umpy.any(numpy.isnan(field)):99 if np.any(np.isnan(field)): 94 100 md = md.checkmessage(options.getfieldvalue('message',\ 95 101 "NaN values found in field '%s'" % fieldname)) … … 97 103 #check Inf 98 104 if options.getfieldvalue('Inf',0): 99 if n umpy.any(numpy.isinf(field)):105 if np.any(np.isinf(field)): 100 106 md = md.checkmessage(options.getfieldvalue('message',\ 101 107 "Inf values found in field '%s'" % fieldname)) … … 124 130 if options.exist('>='): 125 131 lowerbound=options.getfieldvalue('>=') 126 if n umpy.any(field<lowerbound):132 if np.any(field<lowerbound): 127 133 md = md.checkmessage(options.getfieldvalue('message',\ 128 134 "field '%s' should have values above %d" % (fieldname,lowerbound))) 129 135 if options.exist('>'): 130 136 lowerbound=options.getfieldvalue('>') 131 if n umpy.any(field<=lowerbound):137 if np.any(field<=lowerbound): 132 138 md = md.checkmessage(options.getfieldvalue('message',\ 133 139 "field '%s' should have values above %d" % (fieldname,lowerbound))) … … 136 142 if options.exist('<='): 137 143 upperbound=options.getfieldvalue('<=') 138 if n umpy.any(field>upperbound):144 if np.any(field>upperbound): 139 145 md = md.checkmessage(options.getfieldvalue('message',\ 140 146 "field '%s' should have values below %d" % (fieldname,upperbound))) 141 147 if options.exist('<'): 142 148 upperbound=options.getfieldvalue('<') 143 if n umpy.any(field>=upperbound):149 if np.any(field>=upperbound): 144 150 md = md.checkmessage(options.getfieldvalue('message',\ 145 151 "field '%s' should have values below %d" % (fieldname,upperbound))) … … 158 164 #Check forcings (size and times) 159 165 if options.getfieldvalue('timeseries',0): 160 if n umpy.size(field,0)==md.mesh.numberofvertices:161 if n umpy.ndim(field)>1 and not numpy.size(field,1)==1:166 if np.size(field,0)==md.mesh.numberofvertices: 167 if np.ndim(field)>1 and not np.size(field,1)==1: 162 168 md = md.checkmessage(options.getfieldvalue('message',\ 163 169 "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) 164 elif n umpy.size(field,0)==md.mesh.numberofvertices+1 or numpy.size(field,0)==2:165 if not all(field[-1,:]==n umpy.sort(field[-1,:])):170 elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==2: 171 if not all(field[-1,:]==np.sort(field[-1,:])): 166 172 md = md.checkmessage(options.getfieldvalue('message',\ 167 173 "field '%s' columns should be sorted chronologically" % fieldname)) … … 175 181 #Check single value forcings (size and times) 176 182 if options.getfieldvalue('singletimeseries',0): 177 if n umpy.size(field,0)==2:178 if not all(field[-1,:]==n umpy.sort(field[-1,:])):183 if np.size(field,0)==2: 184 if not all(field[-1,:]==np.sort(field[-1,:])): 179 185 md = md.checkmessage(options.getfieldvalue('message',\ 180 186 "field '%s' columns should be sorted chronologically" % fieldname)) -
issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py
r21268 r21303 1 1 from netCDF4 import Dataset, stringtochar 2 import numpy as np2 import numpy as np 3 3 import time 4 4 import collections … … 42 42 typelist=[bool,str,unicode,int,float,complex, 43 43 collections.OrderedDict, 44 n umpy.int64,numpy.ndarray,numpy.float64]44 np.int64,np.ndarray,np.float64] 45 45 groups=dict.keys(md.__dict__) 46 46 for group in groups: … … 121 121 TypeDict = {float:'f8', 122 122 'float64':'f8', 123 n umpy.float64:'f8',123 np.float64:'f8', 124 124 int:'i8', 125 125 'int64':'i8', 126 n umpy.int64:'i8',126 np.int64:'i8', 127 127 str:str, 128 128 dict:str} -
issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py
r20786 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 import model … … 40 40 os.mkdir(filename) 41 41 42 IsEnveloppe=n umpy.where(model.mesh.vertexonbase | model.mesh.vertexonsurface)42 IsEnveloppe=np.where(model.mesh.vertexonbase | model.mesh.vertexonsurface) 43 43 #get the element related variables 44 44 if 'z' in dict.keys(model.mesh.__dict__): 45 points=n umpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z))46 num_of_elt=n umpy.size(numpy.isnan(model.mesh.lowerelements))+numpy.size(numpy.isnan(model.mesh.upperelements))47 low_elt_num=n umpy.size(numpy.isnan(model.mesh.lowerelements))48 top_elt_num=n umpy.size(numpy.isnan(model.mesh.upperelements))45 points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) 46 num_of_elt=np.size(np.isnan(model.mesh.lowerelements))+np.size(np.isnan(model.mesh.upperelements)) 47 low_elt_num=np.size(np.isnan(model.mesh.lowerelements)) 48 top_elt_num=np.size(np.isnan(model.mesh.upperelements)) 49 49 else: 50 points=n umpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x))))51 num_of_elt=n umpy.shape(model.mesh.elements)[0]50 points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) 51 num_of_elt=np.shape(model.mesh.elements)[0] 52 52 53 num_of_points=n umpy.size(points)[0]54 dim=n umpy.size(points)[1]55 point_per_elt=n umpy.shape(model.mesh.elements)[1]53 num_of_points=np.size(points)[0] 54 dim=np.size(points)[1] 55 point_per_elt=np.shape(model.mesh.elements)[1] 56 56 57 57 celltype=5 #triangles … … 67 67 for solution in solnames: 68 68 #looking for multiple time steps 69 if (n umpy.size(res_struct.__dict__[solution])>num_of_timesteps):70 num_of_timesteps=n umpy.size(res_struct.__dict__[solution])69 if (np.size(res_struct.__dict__[solution])>num_of_timesteps): 70 num_of_timesteps=np.size(res_struct.__dict__[solution]) 71 71 num_of_timesteps=int(num_of_timesteps)+1 72 72 else: … … 106 106 for sol in solnames: 107 107 #dealing with results on different timesteps 108 if(n umpy.size(res_struct.__dict__[sol])>timestep):108 if(np.size(res_struct.__dict__[sol])>timestep): 109 109 timestep = step 110 110 else: 111 timestep = n umpy.size(res_struct.__dict__[sol])111 timestep = np.size(res_struct.__dict__[sol]) 112 112 113 113 #getting the fields in the solution 114 if(n umpy.size(res_struct.__dict__[sol])>1):114 if(np.size(res_struct.__dict__[sol])>1): 115 115 fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep-1).__dict__) 116 116 else: … … 118 118 #check which field is a real result and print 119 119 for field in fieldnames: 120 if(n umpy.size(res_struct.__dict__[sol])>1):120 if(np.size(res_struct.__dict__[sol])>1): 121 121 fieldstruct=res_struct.__dict__[sol].__getitem__(timestep-1).__dict__[field] 122 122 else: 123 123 fieldstruct=res_struct.__dict__[sol].__dict__[field] 124 124 125 if ((n umpy.size(fieldstruct))==num_of_points):125 if ((np.size(fieldstruct))==num_of_points): 126 126 fid.write('SCALARS %s float 1 \n' % field) 127 127 fid.write('LOOKUP_TABLE default\n') 128 128 for node in range(0,num_of_points): 129 129 #paraview does not like NaN, replacing 130 if n umpy.isnan(fieldstruct[node]):130 if np.isnan(fieldstruct[node]): 131 131 fid.write('%e\n' % -9999.9999) 132 132 #also checking for verry small value that mess up … … 143 143 for field in othernames: 144 144 # fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe)); 145 if ((n umpy.size(other_struct.__dict__[field]))==num_of_points):145 if ((np.size(other_struct.__dict__[field]))==num_of_points): 146 146 fid.write('SCALARS %s float 1 \n' % field) 147 147 fid.write('LOOKUP_TABLE default\n') 148 148 for node in range(0,num_of_points): 149 149 #paraview does not like NaN, replacing 150 if n umpy.isnan(other_struct.__dict__[field][node]):150 if np.isnan(other_struct.__dict__[field][node]): 151 151 fid.write('%e\n' % -9999.9999) 152 152 #also checking for verry small value that mess up -
issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py
r21130 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 import model … … 42 42 #get the element related variables 43 43 if 'z' in dict.keys(model.mesh.__dict__): 44 points=n umpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z))44 points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) 45 45 dim=3 46 46 else: 47 points=n umpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x))))47 points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) 48 48 dim=2 49 49 50 num_of_points=n umpy.size(model.mesh.x)51 num_of_elt=n umpy.shape(model.mesh.elements)[0]52 point_per_elt=n umpy.shape(model.mesh.elements)[1]50 num_of_points=np.size(model.mesh.x) 51 num_of_elt=np.shape(model.mesh.elements)[0] 52 point_per_elt=np.shape(model.mesh.elements)[1] 53 53 54 54 #Select the type of element function of the number of nodes per elements … … 70 70 for solution in solnames: 71 71 #looking for multiple time steps 72 if (n umpy.size(res_struct.__dict__[solution])>num_of_timesteps):73 num_of_timesteps=n umpy.size(res_struct.__dict__[solution])72 if (np.size(res_struct.__dict__[solution])>num_of_timesteps): 73 num_of_timesteps=np.size(res_struct.__dict__[solution]) 74 74 num_of_timesteps=int(num_of_timesteps) 75 75 else: … … 112 112 for sol in solnames: 113 113 #dealing with results on different timesteps 114 if(n umpy.size(res_struct.__dict__[sol])>timestep):114 if(np.size(res_struct.__dict__[sol])>timestep): 115 115 timestep = step 116 116 else: 117 timestep = n umpy.size(res_struct.__dict__[sol])117 timestep = np.size(res_struct.__dict__[sol]) 118 118 119 119 #getting the fields in the solution 120 if(n umpy.size(res_struct.__dict__[sol])>1):120 if(np.size(res_struct.__dict__[sol])>1): 121 121 fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep).__dict__) 122 122 else: … … 124 124 #check which field is a real result and print 125 125 for field in fieldnames: 126 if(n umpy.size(res_struct.__dict__[sol])>1):126 if(np.size(res_struct.__dict__[sol])>1): 127 127 fieldstruct=res_struct.__dict__[sol].__getitem__(timestep).__dict__[field] 128 128 else: 129 129 fieldstruct=res_struct.__dict__[sol].__dict__[field] 130 130 131 if ((n umpy.size(fieldstruct))==num_of_points):131 if ((np.size(fieldstruct))==num_of_points): 132 132 fid.write('SCALARS %s float 1 \n' % field) 133 133 fid.write('LOOKUP_TABLE default\n') 134 134 for node in range(0,num_of_points): 135 135 #paraview does not like NaN, replacing 136 if n umpy.isnan(fieldstruct[node]):136 if np.isnan(fieldstruct[node]): 137 137 fid.write('%e\n' % -9999.9999) 138 138 #also checking for verry small value that mess up … … 148 148 othernames=(dict.keys(other_struct.__dict__)) 149 149 for field in othernames: 150 if ((n umpy.size(other_struct.__dict__[field]))==num_of_points):150 if ((np.size(other_struct.__dict__[field]))==num_of_points): 151 151 fid.write('SCALARS %s float 1 \n' % field) 152 152 fid.write('LOOKUP_TABLE default\n') 153 153 for node in range(0,num_of_points): 154 154 #paraview does not like NaN, replacing 155 if n umpy.isnan(other_struct.__dict__[field][node]):155 if np.isnan(other_struct.__dict__[field][node]): 156 156 fid.write('%e\n' % -9999.9999) 157 157 #also checking for verry small value that mess up -
issm/trunk-jpl/src/m/contrib/morlighem/bamg/YamsCall.py
r17480 r21303 1 import numpy 1 import numpy as np 2 2 import time 3 3 import subprocess … … 37 37 t1=time.time() 38 38 print "%s" % ' computing metric...' 39 metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,n umpy.empty(0,int))39 metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,np.empty(0,int)) 40 40 t2=time.time() 41 41 print "%s%d%s\n" % (' done (',t2-t1,' seconds)') … … 44 44 t1=time.time() 45 45 print "%s" % ' writing initial mesh files...' 46 n umpy.savetxt('carre0.met',metric)46 np.savetxt('carre0.met',metric) 47 47 48 48 f=open('carre0.mesh','w') … … 66 66 67 67 #Deal with rifts 68 if n umpy.any(not numpy.isnan(md.rifts.riftstruct)):68 if np.any(not np.isnan(md.rifts.riftstruct)): 69 69 70 70 #we have the list of triangles that make up the rift. keep those triangles around during refinement. 71 triangles=n umpy.empty(0,int)71 triangles=np.empty(0,int) 72 72 for riftstruct in md.rifts.riftstruct: 73 triangles=n umpy.concatenate((triangles,riftstruct.segments[:,2]))73 triangles=np.concatenate((triangles,riftstruct.segments[:,2])) 74 74 75 f.write("\n\n%s\n%i\n\n" % ('RequiredTriangles',n umpy.size(triangles)))75 f.write("\n\n%s\n%i\n\n" % ('RequiredTriangles',np.size(triangles))) 76 76 for triangle in triangles: 77 77 f.write("%i\n" % triangle) … … 97 97 t1=time.time() 98 98 print "\n%s" % ' reading final mesh files...' 99 Tria=n umpy.loadtxt('carre1.tria',int)100 Coor=n umpy.loadtxt('carre1.coor',float)99 Tria=np.loadtxt('carre1.tria',int) 100 Coor=np.loadtxt('carre1.coor',float) 101 101 md.mesh.x=Coor[:,0] 102 102 md.mesh.y=Coor[:,1] 103 md.mesh.z=n umpy.zeros((numpy.size(Coor,axis=0),1))103 md.mesh.z=np.zeros((np.size(Coor,axis=0),1)) 104 104 md.mesh.elements=Tria 105 md.mesh.numberofvertices=n umpy.size(Coor,axis=0)106 md.mesh.numberofelements=n umpy.size(Tria,axis=0)105 md.mesh.numberofvertices=np.size(Coor,axis=0) 106 md.mesh.numberofelements=np.size(Tria,axis=0) 107 107 numberofelements2=md.mesh.numberofelements 108 108 t2=time.time() -
issm/trunk-jpl/src/m/coordsystems/gmtmask.py
r21069 r21303 1 1 from MatlabFuncs import * 2 2 from model import * 3 from n umpyimport *3 from np.import * 4 4 from os import getenv, putenv 5 5 import subprocess -
issm/trunk-jpl/src/m/coordsystems/ll2xy.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71): -
issm/trunk-jpl/src/m/coordsystems/xy2ll.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from math import pi 3 3 … … 37 37 raise StandardError('bad usage: type "help(xy2ll)" for details') 38 38 39 # if x,y passed as lists, convert to n umpyarrays40 if type(x) != "n umpy.ndarray":39 # if x,y passed as lists, convert to np.arrays 40 if type(x) != "np.ndarray": 41 41 x=np.array(x) 42 if type(y) != "n umpy.ndarray":42 if type(y) != "np.ndarray": 43 43 y=np.array(y) 44 44 -
issm/trunk-jpl/src/m/exp/expcoarsen.py
r21254 r21303 1 1 import os.path 2 import numpy as np2 import numpy as np 3 3 from collections import OrderedDict 4 4 from expread import expread -
issm/trunk-jpl/src/m/exp/expdisp.py
r21278 r21303 1 1 from expread import expread 2 import numpy as np2 import numpy as np 3 3 from matplotlib.path import Path 4 4 import matplotlib.patches as patches -
issm/trunk-jpl/src/m/exp/expread.py
r17731 r21303 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 from collections import OrderedDict 4 4 import MatlabFuncs as m 5 5 6 6 def expread(filename): 7 7 8 """ 9 8 10 EXPREAD - read a file exp and build a Structure 9 11 … … 23 25 24 26 See also EXPDOC, EXPWRITEASVERTICES 27 25 28 """ 26 27 29 #some checks 28 30 if not os.path.exists(filename): … … 31 33 #initialize number of profile 32 34 contours=[] 33 34 35 #open file 35 36 fid=open(filename,'r') 36 37 37 #loop over the number of profiles 38 38 while True: 39 40 39 #update number of profiles 41 40 contour=OrderedDict() 42 43 41 #Get file name 44 42 A=fid.readline() … … 50 48 if not (len(A) == 2 and m.strcmp(A[0],'##') and m.strncmp(A[1],'Name:',5)): 51 49 break 50 52 51 if len(A[1])>5: 53 52 contour['name']=A[1][5:-1] … … 59 58 if not (len(A) == 2 and m.strcmp(A[0],'##') and m.strncmp(A[1],'Icon:',5)): 60 59 break 61 62 60 #Get Info 63 61 A=fid.readline().split() … … 72 70 #Get Info 73 71 A=fid.readline().split() 74 if not (len(A) == 5 and m.strcmp(A[0],'#') and m.strcmp(A[1],'X') and m.strcmp(A[2],'pos') \75 72 if not (len(A) == 5 and m.strcmp(A[0],'#') and m.strcmp(A[1],'X') and m.strcmp(A[2],'pos') 73 and m.strcmp(A[3],'Y') and m.strcmp(A[4],'pos')): 76 74 break 77 78 75 #Get Coordinates 79 contour['x']=n umpy.empty(contour['nods'])80 contour['y']=n umpy.empty(contour['nods'])76 contour['x']=np.empty(contour['nods']) 77 contour['y']=np.empty(contour['nods']) 81 78 for i in xrange(int(contour['nods'])): 82 79 A=fid.readline().split() … … 93 90 94 91 contours.append(contour) 95 96 92 #close file 97 93 fid.close() 98 99 94 return contours 100 -
issm/trunk-jpl/src/m/exp/expwrite.py
r19459 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def expwrite(contours,filename): … … 21 21 fid=open(filename,'w') 22 22 for x,y in zip(contours['x'],contours['y']): 23 #if n umpy.size(contour['x'])!=numpy.size(contour['y']):23 #if np.size(contour['x'])!=np.size(contour['y']): 24 24 if len(x)!=len(y): 25 25 raise RuntimeError("contours x and y coordinates must be of identical size") … … 36 36 fid.write("%s\n" % '## Icon:0') 37 37 fid.write("%s\n" % '# Points Count Value') 38 #fid.write("%i %f\n" % (n umpy.size(contour['x']),contour['density']))39 fid.write("%i %f\n" % (n umpy.size(x),density))38 #fid.write("%i %f\n" % (np.size(contour['x']),contour['density'])) 39 fid.write("%i %f\n" % (np.size(x),density)) 40 40 fid.write("%s\n" % '# X pos Y pos') 41 41 #for x,y in zip(contour['x'],contour['y']): -
issm/trunk-jpl/src/m/extrusion/DepthAverage.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from project2d import project2d 3 3 … … 46 46 47 47 if vec2d: 48 vector_average=vector_average.reshape(-1, 1)48 vector_average=vector_average.reshape(-1,) 49 49 50 50 return vector_average -
issm/trunk-jpl/src/m/extrusion/project2d.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def project2d(md3d,value,layer): … … 45 45 46 46 if vec2d: 47 projection_value=projection_value.reshape(-1, 1)47 projection_value=projection_value.reshape(-1,) 48 48 49 49 return projection_value -
issm/trunk-jpl/src/m/extrusion/project3d.py
r17725 r21303 1 import numpy 1 import numpy as np 2 2 from pairoptions import pairoptions 3 3 … … 37 37 38 38 vector1d=False 39 if isinstance(vector2d,n umpy.ndarray) and numpy.ndim(vector2d)==1:39 if isinstance(vector2d,np.ndarray) and np.ndim(vector2d)==1: 40 40 vector1d=True 41 vector2d=vector2d.reshape(-1, 1)41 vector2d=vector2d.reshape(-1,) 42 42 43 if isinstance(vector2d,(bool,int,long,float)) or n umpy.size(vector2d)==1:43 if isinstance(vector2d,(bool,int,long,float)) or np.size(vector2d)==1: 44 44 projected_vector=vector2d 45 45 … … 47 47 48 48 #Initialize 3d vector 49 if vector2d.shape[0]==md.mesh.numberofvertices2d: 50 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 51 elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: 52 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 53 projected_vector[-1,:]=vector2d[-1,:] 54 vector2d=vector2d[:-1,:] 49 if np.ndim(vector2d)==1: 50 if vector2d.shape[0]==md.mesh.numberofvertices2d: 51 projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype) 52 elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: 53 projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1))).astype(vector2d.dtype) 54 projected_vector[-1]=vector2d[-1] 55 vector2d=vector2d[:-1] 56 else: 57 raise TypeError("vector length not supported") 58 #Fill in 59 if layer==0: 60 for i in xrange(md.mesh.numberoflayers): 61 projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d)]=vector2d 62 else: 63 projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d)]=vector2d 55 64 else: 56 raise TypeError("vector length not supported") 65 if vector2d.shape[0]==md.mesh.numberofvertices2d: 66 projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices,np.size(vector2d,axis=1)))).astype(vector2d.dtype) 67 elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: 68 projected_vector=(paddingvalue*np.ones((md.mesh.numberofvertices+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) 69 projected_vector[-1,:]=vector2d[-1,:] 70 vector2d=vector2d[:-1,:] 71 else: 72 raise TypeError("vector length not supported") 73 #Fill in 74 if layer==0: 75 for i in xrange(md.mesh.numberoflayers): 76 projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d),:]=vector2d 77 else: 78 projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d 57 79 58 #Fill in59 if layer==0:60 for i in xrange(md.mesh.numberoflayers):61 projected_vector[(i*md.mesh.numberofvertices2d):((i+1)*md.mesh.numberofvertices2d),:]=vector2d62 else:63 projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d64 80 65 81 elif vectype.lower()=='element': 66 82 67 83 #Initialize 3d vector 68 if vector2d.shape[0]==md.mesh.numberofelements2d: 69 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 70 elif vector2d.shape[0]==md.mesh.numberofelements2d+1: 71 projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) 72 projected_vector[-1,:]=vector2d[-1,:] 73 vector2d=vector2d[:-1,:] 84 if np.ndim(vector2d)==1: 85 if vector2d.shape[0]==md.mesh.numberofelements2d: 86 projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements))).astype(vector2d.dtype) 87 elif vector2d.shape[0]==md.mesh.numberofelements2d+1: 88 projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1))).astype(vector2d.dtype) 89 projected_vector[-1]=vector2d[-1] 90 vector2d=vector2d[:-1] 91 else: 92 raise TypeError("vector length not supported") 93 #Fill in 94 if layer==0: 95 for i in xrange(md.mesh.numberoflayers-1): 96 projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d)]=vector2d 97 else: 98 projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d)]=vector2d 74 99 else: 75 raise TypeError("vector length not supported") 76 77 #Fill in 78 if layer==0: 79 for i in xrange(md.mesh.numberoflayers-1): 80 projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d),:]=vector2d 81 else: 82 projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d),:]=vector2d 100 if vector2d.shape[0]==md.mesh.numberofelements2d: 101 projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements, np.size(vector2d,axis=1)))).astype(vector2d.dtype) 102 elif vector2d.shape[0]==md.mesh.numberofelements2d+1: 103 projected_vector=(paddingvalue*np.ones((md.mesh.numberofelements+1,np.size(vector2d,axis=1)))).astype(vector2d.dtype) 104 projected_vector[-1,:]=vector2d[-1,:] 105 vector2d=vector2d[:-1,:] 106 else: 107 raise TypeError("vector length not supported") 108 #Fill in 109 if layer==0: 110 for i in xrange(md.mesh.numberoflayers-1): 111 projected_vector[(i*md.mesh.numberofelements2d):((i+1)*md.mesh.numberofelements2d),:]=vector2d 112 else: 113 projected_vector[((layer-1)*md.mesh.numberofelements2d):(layer*md.mesh.numberofelements2d),:]=vector2d 83 114 84 115 else: -
issm/trunk-jpl/src/m/geometry/FlagElements.py
r20910 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 #from basinzoom import basinzoon … … 24 24 if isinstance(region,(str,unicode)): 25 25 if not region: 26 flag=n umpy.zeros(md.mesh.numberofelements,bool)26 flag=np.zeros(md.mesh.numberofelements,bool) 27 27 invert=0 28 28 elif m.strcmpi(region,'all'): 29 flag=n umpy.ones(md.mesh.numberofelements,bool)29 flag=np.ones(md.mesh.numberofelements,bool) 30 30 invert=0 31 31 else: … … 44 44 xlim,ylim=basinzoom('basin',region) 45 45 flag_nodes=p.logical_and_n(md.mesh.x<xlim[1],md.mesh.x>xlim[0],md.mesh.y<ylim[1],md.mesh.y>ylim[0]) 46 flag=n umpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool)46 flag=np.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool) 47 47 else: 48 48 #ok, flag elements … … 51 51 52 52 if invert: 53 flag=n umpy.logical_not(flag)53 flag=np.logical_not(flag) 54 54 55 elif isinstance(region,n umpy.ndarray) or isinstance(region,bool):56 if n umpy.size(region,0)==md.mesh.numberofelements:55 elif isinstance(region,np.ndarray) or isinstance(region,bool): 56 if np.size(region,0)==md.mesh.numberofelements: 57 57 flag=region 58 elif n umpy.size(region,0)==md.mesh.numberofvertices:59 flag=(n umpy.sum(region[md.mesh.elements-1]>0,axis=1)==numpy.size(md.mesh.elements,1))58 elif np.size(region,0)==md.mesh.numberofvertices: 59 flag=(np.sum(region[md.mesh.elements-1]>0,axis=1)==np.size(md.mesh.elements,1)) 60 60 else: 61 61 raise TypeError("Flaglist for region must be of same size as number of elements in model.") -
issm/trunk-jpl/src/m/geometry/GetAreas.py
r13984 r21303 1 import numpy 1 import numpy as np 2 2 3 def GetAreas(index,x,y,z=n umpy.array([])):3 def GetAreas(index,x,y,z=np.array([])): 4 4 """ 5 5 GETAREAS - compute areas or volumes of elements … … 18 18 19 19 #get number of elements and number of nodes 20 nels=n umpy.size(index,axis=0)21 nods=n umpy.size(x)20 nels=np.size(index,axis=0) 21 nods=np.size(x) 22 22 23 23 #some checks 24 if n umpy.size(y)!=nods or (z and numpy.size(z)!=nods):24 if np.size(y)!=nods or (z and np.size(z)!=nods): 25 25 raise TypeError("GetAreas error message: x,y and z do not have the same length.") 26 if n umpy.max(index)>nods:26 if np.max(index)>nods: 27 27 raise TypeError("GetAreas error message: index should not have values above %d." % nods) 28 if (not z and n umpy.size(index,axis=1)!=3):28 if (not z and np.size(index,axis=1)!=3): 29 29 raise TypeError("GetAreas error message: index should have 3 columns for 2d meshes.") 30 if (z and n umpy.size(index,axis=1)!=6):30 if (z and np.size(index,axis=1)!=6): 31 31 raise TypeError("GetAreas error message: index should have 6 columns for 3d meshes.") 32 32 33 33 #initialization 34 areas=n umpy.zeros(nels)34 areas=np.zeros(nels) 35 35 x1=x[index[:,0]-1] 36 36 x2=x[index[:,1]-1] … … 46 46 else: 47 47 #V=area(triangle)*1/3(z1+z2+z3) 48 thickness=n umpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1])48 thickness=np.mean(z[index[:,3:6]-1])-np.mean(z[index[:,0:3]-1]) 49 49 areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness 50 50 -
issm/trunk-jpl/src/m/geometry/SegIntersect.py
r13503 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def SegIntersect(seg1,seg2): … … 24 24 yD=seg2[1,1] 25 25 26 O2A=n umpy.array([xA,yA])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.])27 O2B=n umpy.array([xB,yB])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.])28 O1C=n umpy.array([xC,yC])-numpy.array([xA/2.+xB/2.,yB/2.+yA/2.])29 O1D=n umpy.array([xD,yD])-numpy.array([xA/2.+xB/2.,yB/2.+yA/2.])26 O2A=np.array([xA,yA])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) 27 O2B=np.array([xB,yB])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) 28 O1C=np.array([xC,yC])-np.array([xA/2.+xB/2.,yB/2.+yA/2.]) 29 O1D=np.array([xD,yD])-np.array([xA/2.+xB/2.,yB/2.+yA/2.]) 30 30 31 n1=n umpy.array([yA-yB,xB-xA]) #normal vector to segA32 n2=n umpy.array([yC-yD,xD-xC]) #normal vector to segB31 n1=np.array([yA-yB,xB-xA]) #normal vector to segA 32 n2=np.array([yC-yD,xD-xC]) #normal vector to segB 33 33 34 test1=n umpy.dot(n2,O2A)35 test2=n umpy.dot(n2,O2B)34 test1=np.dot(n2,O2A) 35 test2=np.dot(n2,O2B) 36 36 37 37 if test1*test2>0: … … 39 39 return bval 40 40 41 test3=n umpy.dot(n1,O1C)42 test4=n umpy.dot(n1,O1D)41 test3=np.dot(n1,O1C) 42 test4=np.dot(n1,O1D) 43 43 44 44 if test3*test4>0: … … 47 47 48 48 #if colinear 49 if test1*test2==0 and test3*test4==0 and n umpy.linalg.det(numpy.hstack((n1.reshape((-1,1)),n2.reshape(-1,1))))==0:49 if test1*test2==0 and test3*test4==0 and np.linalg.det(np.hstack((n1.reshape((-1,)),n2.reshape(-1,))))==0: 50 50 51 51 #projection on the axis O1O2 52 O2O1=n umpy.array([xA/2.+xB/2.,yB/2.+yA/2.])-numpy.array([xD/2.+xC/2.,yD/2.+yC/2.])53 O1A=n umpy.dot(O2O1,(O2A-O2O1))54 O1B=n umpy.dot(O2O1,(O2B-O2O1))55 O1C=n umpy.dot(O2O1,O1C)56 O1D=n umpy.dot(O2O1,O1D)52 O2O1=np.array([xA/2.+xB/2.,yB/2.+yA/2.])-np.array([xD/2.+xC/2.,yD/2.+yC/2.]) 53 O1A=np.dot(O2O1,(O2A-O2O1)) 54 O1B=np.dot(O2O1,(O2B-O2O1)) 55 O1C=np.dot(O2O1,O1C) 56 O1D=np.dot(O2O1,O1D) 57 57 58 58 #test if one point is included in the other segment (->bval=1) -
issm/trunk-jpl/src/m/geometry/slope.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff 3 3 -
issm/trunk-jpl/src/m/interp/SectionValues.py
r21254 r21303 1 1 import os 2 2 from expread import expread 3 import numpy as np3 import numpy as np 4 4 from project2d import project2d 5 5 #from InterpFromMesh2d import InterpFromMesh2d -
issm/trunk-jpl/src/m/interp/averaging.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from GetAreas import GetAreas 3 3 import MatlabFuncs as m … … 82 82 average_node = csc_matrix(average_node) 83 83 else: 84 average_node=csc_matrix(data.reshape(-1, 1))84 average_node=csc_matrix(data.reshape(-1,)) 85 85 86 86 #loop over iteration -
issm/trunk-jpl/src/m/interp/holefiller.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from scipy.spatial import cKDTree 3 3 -
issm/trunk-jpl/src/m/interp/interp.py
r21254 r21303 1 1 # module for inperpolating/smoothing data 2 import numpy as np2 import numpy as np 3 3 from scipy.interpolate import CloughTocher2DInterpolator, Rbf 4 4 from scipy.spatial import cKDTree -
issm/trunk-jpl/src/m/io/loadvars.py
r21268 r21303 1 1 import shelve 2 2 import os.path 3 import numpy as np3 import numpy as np 4 4 from netCDF4 import Dataset 5 5 from netCDF4 import chartostring … … 7 7 from os import path 8 8 from whichdb import whichdb 9 from collections import OrderedDict 9 10 from model import * 10 11 … … 130 131 Tree[t].__dict__[str(var)]=varval[0] 131 132 else: 132 if varval[0]==u'':133 if str(varval[0])=='': 133 134 Tree.__dict__[str(var)]=[] 134 135 elif varval[0]=='True': … … 167 168 #dealling with dict 168 169 if varval.dtype==str: 169 Tree.__dict__[str(var)]= dict(zip(varval[:,0], varval[:,1]))170 Tree.__dict__[str(var)]=OrderedDict(zip(varval[:,0], varval[:,1])) 170 171 else: 171 172 if type(Tree)==list: … … 178 179 else: 179 180 Tree.__dict__[str(var)]=varval[:,:] 181 try: 182 print('treated {}.{} with type {}'.format(mod,var,type(Tree.__dict__[str(var)]))) 183 except AttributeError: 184 print('treated {}.{} with type {}'.format(mod,var,type(Tree[0].__dict__[str(var)]))) 180 185 elif vardim==3: 181 186 if type(Tree)==list: -
issm/trunk-jpl/src/m/materials/TMeltingPoint.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def TMeltingPoint(reftemp,pressure): -
issm/trunk-jpl/src/m/materials/cuffey.py
r20628 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def cuffey(temperature): … … 13 13 """ 14 14 15 if n umpy.any(temperature<0.):15 if np.any(temperature<0.): 16 16 raise RuntimeError("input temperature should be in Kelvin (positive)") 17 17 18 18 T = temperature-273.15 19 19 20 rigidity=n umpy.zeros_like(T)21 pos=n umpy.nonzero(T<=-45)20 rigidity=np.zeros_like(T) 21 pos=np.nonzero(T<=-45) 22 22 rigidity[pos]=10**8*(-0.000396645116301*(T[pos]+50)**3+ 0.013345579471334*(T[pos]+50)**2 -0.356868703259105*(T[pos]+50)+7.272363035371383) 23 pos=n umpy.nonzero(numpy.logical_and(-45<=T,T<-40))23 pos=np.nonzero(np.logical_and(-45<=T,T<-40)) 24 24 rigidity[pos]=10**8*(-0.000396645116301*(T[pos]+45)**3+ 0.007395902726819*(T[pos]+45)**2 -0.253161292268336*(T[pos]+45)+5.772078366321591) 25 pos=n umpy.nonzero(numpy.logical_and(-40<=T,T<-35))25 pos=np.nonzero(np.logical_and(-40<=T,T<-35)) 26 26 rigidity[pos]=10**8*(0.000408322072669*(T[pos]+40)**3+ 0.001446225982305*(T[pos]+40)**2 -0.208950648722716*(T[pos]+40)+4.641588833612773) 27 pos=n umpy.nonzero(numpy.logical_and(-35<=T,T<-30))27 pos=np.nonzero(np.logical_and(-35<=T,T<-30)) 28 28 rigidity[pos]=10**8*(-0.000423888728124*(T[pos]+35)**3+ 0.007571057072334*(T[pos]+35)**2 -0.163864233449525*(T[pos]+35)+3.684031498640382) 29 pos=n umpy.nonzero(numpy.logical_and(-30<=T,T<-25))29 pos=np.nonzero(np.logical_and(-30<=T,T<-25)) 30 30 rigidity[pos]=10**8*(0.000147154327025*(T[pos]+30)**3+ 0.001212726150476*(T[pos]+30)**2 -0.119945317335478*(T[pos]+30)+3.001000667185614) 31 pos=n umpy.nonzero(numpy.logical_and(-25<=T,T<-20))31 pos=np.nonzero(np.logical_and(-25<=T,T<-20)) 32 32 rigidity[pos]=10**8*(-0.000193435838672*(T[pos]+25)**3+ 0.003420041055847*(T[pos]+25)**2 -0.096781481303861*(T[pos]+25)+2.449986525148220) 33 pos=n umpy.nonzero(numpy.logical_and(-20<=T,T<-15))33 pos=np.nonzero(np.logical_and(-20<=T,T<-15)) 34 34 rigidity[pos]=10**8*(0.000219771255067*(T[pos]+20)**3+ 0.000518503475772*(T[pos]+20)**2 -0.077088758645767*(T[pos]+20)+2.027400665191131) 35 pos=n umpy.nonzero(numpy.logical_and(-15<=T,T<-10))35 pos=np.nonzero(np.logical_and(-15<=T,T<-10)) 36 36 rigidity[pos]=10**8*(-0.000653438900191*(T[pos]+15)**3+ 0.003815072301777*(T[pos]+15)**2 -0.055420879758021*(T[pos]+15)+1.682390865739973) 37 pos=n umpy.nonzero(numpy.logical_and(-10<=T,T<-5))37 pos=np.nonzero(np.logical_and(-10<=T,T<-5)) 38 38 rigidity[pos]=10**8*(0.000692439419762*(T[pos]+10)**3 -0.005986511201093 *(T[pos]+10)**2 -0.066278074254598*(T[pos]+10)+1.418983411970382) 39 pos=n umpy.nonzero(numpy.logical_and(-5<=T,T<-2))39 pos=np.nonzero(np.logical_and(-5<=T,T<-2)) 40 40 rigidity[pos]=10**8*(-0.000132282004110*(T[pos]+5)**3 +0.004400080095332*(T[pos]+5)**2 -0.074210229783403*(T[pos]+5)+ 1.024485188140279) 41 pos=n umpy.nonzero(-2<=T)41 pos=np.nonzero(-2<=T) 42 42 rigidity[pos]=10**8*(-0.000132282004110*(T[pos]+2)**3 +0.003209542058346*(T[pos]+2)**2 -0.051381363322371*(T[pos]+2)+ 0.837883605537096) 43 43 44 44 #Now make sure that rigidity is positive 45 pos=n umpy.nonzero(rigidity<0)45 pos=np.nonzero(rigidity<0) 46 46 rigidity[pos]=1**6 47 47 -
issm/trunk-jpl/src/m/materials/cuffeytemperate.py
r20628 r21303 1 import numpy 1 import numpy as np 2 2 import cuffey 3 3 … … 15 15 """ 16 16 17 if n umpy.any(temperature<0.):17 if np.any(temperature<0.): 18 18 raise RuntimeError("input temperature should be in Kelvin (positive)") 19 19 20 if (n umpy.any(temperature.shape~=waterfraction.shape)),20 if (np.any(temperature.shape~=waterfraction.shape)), 21 21 error('input temperature and waterfraction should have same size!'); 22 22 end 23 if n umpy.any(waterfraction<0 | waterfraction>1)23 if np.any(waterfraction<0 | waterfraction>1) 24 24 error('input waterfraction should be between 0 and 1'); 25 25 end 26 26 27 rigidity=n umpy.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));27 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)); 28 28 29 29 return rigidity -
issm/trunk-jpl/src/m/materials/paterson.py
r17864 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def paterson(temperature): … … 12 12 """ 13 13 14 if n umpy.any(temperature<0.):14 if np.any(temperature<0.): 15 15 raise RuntimeError("input temperature should be in Kelvin (positive)") 16 16 17 if n umpy.ndim(temperature)==2:17 if np.ndim(temperature)==2: 18 18 #T = temperature.reshape(-1,)-273.15 19 19 T = temperature.flatten()-273.15 20 20 elif isinstance(temperature,float) or isinstance(temperature,int): 21 T = n umpy.array([temperature])-273.1521 T = np.array([temperature])-273.15 22 22 else: 23 23 T = temperature-273.15 … … 36 36 # rigidity=fittedmodel(temperature); 37 37 38 rigidity=n umpy.zeros_like(T)39 pos1=n umpy.nonzero(T<=-45)[0]38 rigidity=np.zeros_like(T) 39 pos1=np.nonzero(T<=-45)[0] 40 40 if len(pos1): 41 41 rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) 42 pos2=n umpy.nonzero(numpy.logical_and(-45<=T,T<-40))[0]42 pos2=np.nonzero(np.logical_and(-45<=T,T<-40))[0] 43 43 if len(pos2): 44 44 rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) 45 pos3=n umpy.nonzero(numpy.logical_and(-40<=T,T<-35))[0]45 pos3=np.nonzero(np.logical_and(-40<=T,T<-35))[0] 46 46 if len(pos3): 47 47 rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) 48 pos4=n umpy.nonzero(numpy.logical_and(-35<=T,T<-30))[0]48 pos4=np.nonzero(np.logical_and(-35<=T,T<-30))[0] 49 49 if len(pos4): 50 50 rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) 51 pos5=n umpy.nonzero(numpy.logical_and(-30<=T,T<-25))[0]51 pos5=np.nonzero(np.logical_and(-30<=T,T<-25))[0] 52 52 if len(pos5): 53 53 rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) 54 pos6=n umpy.nonzero(numpy.logical_and(-25<=T,T<-20))[0]54 pos6=np.nonzero(np.logical_and(-25<=T,T<-20))[0] 55 55 if len(pos6): 56 56 rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) 57 pos7=n umpy.nonzero(numpy.logical_and(-20<=T,T<-15))[0]57 pos7=np.nonzero(np.logical_and(-20<=T,T<-15))[0] 58 58 if len(pos7): 59 59 rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) 60 pos8=n umpy.nonzero(numpy.logical_and(-15<=T,T<-10))[0]60 pos8=np.nonzero(np.logical_and(-15<=T,T<-10))[0] 61 61 if len(pos8): 62 62 rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) 63 pos9=n umpy.nonzero(numpy.logical_and(-10<=T,T<-5))[0]63 pos9=np.nonzero(np.logical_and(-10<=T,T<-5))[0] 64 64 if len(pos9): 65 65 rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) 66 pos10=n umpy.nonzero(numpy.logical_and(-5<=T,T<-2))[0]66 pos10=np.nonzero(np.logical_and(-5<=T,T<-2))[0] 67 67 if len(pos10): 68 68 rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) 69 pos11=n umpy.nonzero(-2<=T)[0]69 pos11=np.nonzero(-2<=T)[0] 70 70 if len(pos11): 71 71 rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) 72 72 73 73 #Now make sure that rigidity is positive 74 pos=n umpy.nonzero(rigidity<0)[0]74 pos=np.nonzero(rigidity<0)[0] 75 75 if len(pos): 76 76 rigidity[pos]=1.e6 -
issm/trunk-jpl/src/m/mech/analyticaldamage.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from averaging import averaging 3 3 #from plotmodel import plotmodel -
issm/trunk-jpl/src/m/mech/backstressfrominversion.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from averaging import averaging 3 3 from thomasparams import thomasparams -
issm/trunk-jpl/src/m/mech/calcbackstress.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from averaging import averaging 3 3 from thomasparams import thomasparams -
issm/trunk-jpl/src/m/mech/damagefrominversion.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def damagefrominversion(md): -
issm/trunk-jpl/src/m/mech/mechanicalproperties.py
r21275 r21303 1 import numpy as np1 import numpy as np 2 2 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff 3 3 from results import results -
issm/trunk-jpl/src/m/mech/robintemperature.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from scipy.special import erf 3 3 -
issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def steadystateiceshelftemp(md,surfacetemp,basaltemp): … … 18 18 temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp) 19 19 """ 20 21 20 if len(md.geometry.thickness)!=md.mesh.numberofvertices: 22 21 raise ValueError('steadystateiceshelftemp error message: thickness should have a length of ' + md.mesh.numberofvertices) -
issm/trunk-jpl/src/m/mech/thomasparams.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from averaging import averaging 3 3 -
issm/trunk-jpl/src/m/mesh/ComputeHessian.py
r17488 r21303 1 import numpy 1 import numpy as np 2 2 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff 3 3 from GetAreas import GetAreas … … 20 20 21 21 #some variables 22 numberofnodes=n umpy.size(x)23 numberofelements=n umpy.size(index,axis=0)22 numberofnodes=np.size(x) 23 numberofelements=np.size(index,axis=0) 24 24 25 25 #some checks 26 if n umpy.size(field)!=numberofnodes and numpy.size(field)!=numberofelements:26 if np.size(field)!=numberofnodes and np.size(field)!=numberofelements: 27 27 raise TypeError("ComputeHessian error message: the given field size not supported yet") 28 28 if not m.strcmpi(type,'node') and not m.strcmpi(type,'element'): … … 38 38 39 39 #compute weights that hold the volume of all the element holding the node i 40 weights=m.sparse(line,n umpy.ones((linesize,1)),numpy.tile(areas.reshape(-1,1),(3,1)),numberofnodes,1)40 weights=m.sparse(line,np.ones((linesize,1)),np.tile(areas.reshape(-1,),(3,1)),numberofnodes,1) 41 41 42 42 #compute field on nodes if on elements 43 if n umpy.size(field,axis=0)==numberofelements:44 field=m.sparse(line,n umpy.ones((linesize,1)),numpy.tile(areas*field,(3,1)),numberofnodes,1)/weights43 if np.size(field,axis=0)==numberofelements: 44 field=m.sparse(line,np.ones((linesize,1)),np.tile(areas*field,(3,1)),numberofnodes,1)/weights 45 45 46 46 #Compute gradient for each element 47 grad_elx=n umpy.sum(field[index-1,0]*alpha,axis=1)48 grad_ely=n umpy.sum(field[index-1,0]*beta,axis=1)47 grad_elx=np.sum(field[index-1,0]*alpha,axis=1) 48 grad_ely=np.sum(field[index-1,0]*beta,axis=1) 49 49 50 50 #Compute gradient for each node (average of the elements around) 51 gradx=m.sparse(line,n umpy.ones((linesize,1)),numpy.tile((areas*grad_elx).reshape(-1,1),(3,1)),numberofnodes,1)52 grady=m.sparse(line,n umpy.ones((linesize,1)),numpy.tile((areas*grad_ely).reshape(-1,1),(3,1)),numberofnodes,1)51 gradx=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_elx).reshape(-1,),(3,1)),numberofnodes,1) 52 grady=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_ely).reshape(-1,),(3,1)),numberofnodes,1) 53 53 gradx=gradx/weights 54 54 grady=grady/weights 55 55 56 56 #Compute hessian for each element 57 hessian=n umpy.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)))57 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 58 58 59 59 if m.strcmpi(type,'node'): 60 60 #Compute Hessian on the nodes (average of the elements around) 61 hessian=n umpy.hstack((m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,0]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \62 m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,1]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \63 m.sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,2]).reshape(-1,1),(3,1)),numberofnodes,1)/weights ))61 hessian=np.hstack((m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,0]).reshape(-1,),(3,1)),numberofnodes,1)/weights, 62 m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,1]).reshape(-1,),(3,1)),numberofnodes,1)/weights, 63 m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,2]).reshape(-1,),(3,1)),numberofnodes,1)/weights )) 64 64 65 65 return hessian -
issm/trunk-jpl/src/m/mesh/ComputeMetric.py
r15749 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos): … … 17 17 b=hessian[:,1] 18 18 d=hessian[:,2] 19 lambda1=0.5*((a+d)+n umpy.sqrt(4.*b**2+(a-d)**2))20 lambda2=0.5*((a+d)-n umpy.sqrt(4.*b**2+(a-d)**2))21 pos1=n umpy.nonzero(lambda1==0.)[0]22 pos2=n umpy.nonzero(lambda2==0.)[0]23 pos3=n umpy.nonzero(numpy.logical_and(b==0.,lambda1==lambda2))[0]19 lambda1=0.5*((a+d)+np.sqrt(4.*b**2+(a-d)**2)) 20 lambda2=0.5*((a+d)-np.sqrt(4.*b**2+(a-d)**2)) 21 pos1=np.nonzero(lambda1==0.)[0] 22 pos2=np.nonzero(lambda2==0.)[0] 23 pos3=np.nonzero(np.logical_and(b==0.,lambda1==lambda2))[0] 24 24 25 25 #Modify the eigen values to control the shape of the elements 26 lambda1=n umpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2)27 lambda2=n umpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2)26 lambda1=np.minimum(np.maximum(np.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) 27 lambda2=np.minimum(np.maximum(np.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) 28 28 29 29 #compute eigen vectors 30 norm1=n umpy.sqrt(8.*b**2+2.*(d-a)**2+2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2))30 norm1=np.sqrt(8.*b**2+2.*(d-a)**2+2.*(d-a)*np.sqrt((a-d)**2+4.*b**2)) 31 31 v1x=2.*b/norm1 32 v1y=((d-a)+n umpy.sqrt((a-d)**2+4.*b**2))/norm133 norm2=n umpy.sqrt(8.*b**2+2.*(d-a)**2-2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2))32 v1y=((d-a)+np.sqrt((a-d)**2+4.*b**2))/norm1 33 norm2=np.sqrt(8.*b**2+2.*(d-a)**2-2.*(d-a)*np.sqrt((a-d)**2+4.*b**2)) 34 34 v2x=2.*b/norm2 35 v2y=((d-a)-n umpy.sqrt((a-d)**2+4.*b**2))/norm235 v2y=((d-a)-np.sqrt((a-d)**2+4.*b**2))/norm2 36 36 37 37 v1x[pos3]=1. … … 41 41 42 42 #Compute new metric (for each node M=V*Lambda*V^-1) 43 metric=n umpy.hstack((((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v2y*v1x-lambda2*v1y*v2x)).reshape(-1,1), \44 ((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v1y*v2y-lambda2*v1y*v2y)).reshape(-1,1), \45 ((v1x*v2y-v1y*v2x)**(-1)*(-lambda1*v2x*v1y+lambda2*v1x*v2y)).reshape(-1,1)))43 metric=np.vstack((((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v2y*v1x-lambda2*v1y*v2x)).reshape(-1,), 44 ((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v1y*v2y-lambda2*v1y*v2y)).reshape(-1,), 45 ((v1x*v2y-v1y*v2x)**(-1)*(-lambda1*v2x*v1y+lambda2*v1x*v2y)).reshape(-1,))).T 46 46 47 47 #some corrections for 0 eigen values 48 metric[pos1,:]=n umpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos1),1))49 metric[pos2,:]=n umpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos2),1))48 metric[pos1,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos1),1)) 49 metric[pos2,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos2),1)) 50 50 51 51 #take care of water elements 52 metric[pos ,:]=n umpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos ),1))52 metric[pos ,:]=np.tile(np.array([[1./hmax**2,0.,1./hmax**2]]),(np.size(pos ),1)) 53 53 54 54 #take care of NaNs if any (use Numpy eig in a loop) 55 pos=n umpy.nonzero(numpy.isnan(metric))[0]56 if n umpy.size(pos):57 print(" %i NaN found in the metric. Use Numpy routine..." % n umpy.size(pos))55 pos=np.nonzero(np.isnan(metric))[0] 56 if np.size(pos): 57 print(" %i NaN found in the metric. Use Numpy routine..." % np.size(pos)) 58 58 for posi in pos: 59 H=n umpy.array([[hessian[posi,0],hessian[posi,1]],[hessian[posi,1],hessian[posi,2]]])60 [v,u]=n umpy.linalg.eig(H)61 v=n umpy.diag(v)59 H=np.array([[hessian[posi,0],hessian[posi,1]],[hessian[posi,1],hessian[posi,2]]]) 60 [v,u]=np.linalg.eig(H) 61 v=np.diag(v) 62 62 lambda1=v[0,0] 63 63 lambda2=v[1,1] 64 v[0,0]=n umpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2)65 v[1,1]=n umpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2)64 v[0,0]=np.minimum(np.maximum(np.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) 65 v[1,1]=np.minimum(np.maximum(np.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) 66 66 67 metricTria=n umpy.dot(numpy.dot(u,v),numpy.linalg.inv(u))68 metric[posi,:]=n umpy.array([metricTria[0,0],metricTria[0,1],metricTria[1,1]])67 metricTria=np.dot(np.dot(u,v),np.linalg.inv(u)) 68 metric[posi,:]=np.array([metricTria[0,0],metricTria[0,1],metricTria[1,1]]) 69 69 70 if n umpy.any(numpy.isnan(metric)):70 if np.any(np.isnan(metric)): 71 71 raise RunTimeError("ComputeMetric error message: NaN in the metric despite our efforts...") 72 72 -
issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py
r17497 r21303 1 import numpy 1 import numpy as np 2 2 import PythonFuncs as p 3 3 … … 12 12 """ 13 13 14 edgeelements=n umpy.nonzero(\15 p.logical_or_n(n umpy.logical_and(elements[:,0]==A,elements[:,1]==B), \16 n umpy.logical_and(elements[:,0]==A,elements[:,2]==B), \17 n umpy.logical_and(elements[:,1]==A,elements[:,2]==B), \18 n umpy.logical_and(elements[:,1]==A,elements[:,0]==B), \19 n umpy.logical_and(elements[:,2]==A,elements[:,0]==B), \20 n umpy.logical_and(elements[:,2]==A,elements[:,1]==B), \14 edgeelements=np.nonzero(\ 15 p.logical_or_n(np.logical_and(elements[:,0]==A,elements[:,1]==B), \ 16 np.logical_and(elements[:,0]==A,elements[:,2]==B), \ 17 np.logical_and(elements[:,1]==A,elements[:,2]==B), \ 18 np.logical_and(elements[:,1]==A,elements[:,0]==B), \ 19 np.logical_and(elements[:,2]==A,elements[:,0]==B), \ 20 np.logical_and(elements[:,2]==A,elements[:,1]==B), \ 21 21 ))[0]+1 22 22 -
issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py
r13984 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def GetNodalFunctionsCoeff(index,x,y): … … 23 23 24 24 #get nels and nods 25 nels=n umpy.size(index,axis=0)26 nods=n umpy.size(x)25 nels=np.size(index,axis=0) 26 nods=np.size(x) 27 27 28 28 #some checks 29 if n umpy.size(y)!=nods:29 if np.size(y)!=nods: 30 30 raise TypeError("GetNodalFunctionsCoeff error message: x and y do not have the same length.") 31 if n umpy.max(index)>nods:31 if np.max(index)>nods: 32 32 raise TypeError("GetNodalFunctionsCoeff error message: index should not have values above %d." % nods) 33 if n umpy.size(index,axis=1)!=3:33 if np.size(index,axis=1)!=3: 34 34 raise TypeError("GetNodalFunctionsCoeff error message: only 2d meshes supported. index should have 3 columns.") 35 35 36 36 #initialize output 37 alpha=n umpy.zeros((nels,3))38 beta=n umpy.zeros((nels,3))37 alpha=np.zeros((nels,3)) 38 beta=np.zeros((nels,3)) 39 39 40 40 #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma … … 48 48 49 49 #get alpha and beta 50 alpha=n umpy.hstack(((invdet*(y2-y3)).reshape(-1,1),(invdet*(y3-y1)).reshape(-1,1),(invdet*(y1-y2)).reshape(-1,1)))51 beta =n umpy.hstack(((invdet*(x3-x2)).reshape(-1,1),(invdet*(x1-x3)).reshape(-1,1),(invdet*(x2-x1)).reshape(-1,1)))50 alpha=np.vstack(((invdet*(y2-y3)).reshape(-1,),(invdet*(y3-y1)).reshape(-1,),(invdet*(y1-y2)).reshape(-1,))).T 51 beta =np.vstack(((invdet*(x3-x2)).reshape(-1,),(invdet*(x1-x3)).reshape(-1,),(invdet*(x2-x1)).reshape(-1,))).T 52 52 53 53 #get gamma if requested 54 gamma=n umpy.zeros((nels,3))55 gamma=n umpy.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)))54 gamma=np.zeros((nels,3)) 55 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 56 56 57 57 return alpha,beta,gamma -
issm/trunk-jpl/src/m/mesh/bamg.py
r20910 r21303 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 from mesh2d import mesh2d 4 4 from collections import OrderedDict … … 94 94 if i: 95 95 flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0] 96 if n umpy.any(numpy.logical_not(flags)):96 if np.any(np.logical_not(flags)): 97 97 raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") 98 98 99 99 #Add all points to bamg_geometry 100 100 nods=domaini['nods']-1 #the domain are closed 0=end 101 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,numpy.hstack((domaini['x'][0:nods].reshape(-1,1),domaini['y'][0:nods].reshape(-1,1),numpy.ones((nods,1))))))102 bamg_geometry.Edges =n umpy.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))))))101 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T)) 102 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)) 103 103 if i: 104 bamg_geometry.SubDomains=numpy.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) 104 bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) 105 106 # 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)))))) 107 # 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,)))))) 108 # if i: 109 # bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) 105 110 106 111 #update counter … … 120 125 #detect whether all points of the rift are inside the domain 121 126 flags=ContourToNodes(rifti['x'],rifti['y'],domain[0],0)[0] 122 if n umpy.all(numpy.logical_not(flags)):127 if np.all(np.logical_not(flags)): 123 128 raise RuntimeError("one rift has all its points outside of the domain outline") 124 129 125 elif n umpy.any(numpy.logical_not(flags)):130 elif np.any(np.logical_not(flags)): 126 131 #We LOTS of work to do 127 132 print "Rift tip outside of or on the domain has been detected and is being processed..." 128 133 129 134 #check that only one point is outside (for now) 130 if n umpy.sum(numpy.logical_not(flags).astype(int))!=1:135 if np.sum(np.logical_not(flags).astype(int))!=1: 131 136 raise RuntimeError("bamg error message: only one point outside of the domain is supported yet") 132 137 … … 136 141 pass 137 142 elif not flags[-1]: 138 rifti['x']=n umpy.flipud(rifti['x'])139 rifti['y']=n umpy.flipud(rifti['y'])143 rifti['x']=np.flipud(rifti['x']) 144 rifti['y']=np.flipud(rifti['y']) 140 145 else: 141 146 raise RuntimeError("bamg error message: only a rift tip can be outside of the domain") … … 146 151 x2=rifti['x'][1] 147 152 y2=rifti['y'][1] 148 for j in xrange(0,n umpy.size(domain[0]['x'])-1):149 if SegIntersect(n umpy.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]]])):153 for j in xrange(0,np.size(domain[0]['x'])-1): 154 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]]])): 150 155 151 156 #Get position of the two nodes of the edge in domain … … 161 166 # 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]); 162 167 # 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]); 163 x=n umpy.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]]))164 y=n umpy.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]]))168 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]])) 169 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]])) 165 170 166 171 segdis= sqrt((x4-x3)**2+(y4-y3)**2) 167 tipdis=n umpy.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)])168 169 if n umpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0):172 tipdis=np.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)]) 173 174 if np.min(tipdis)/segdis < options.getfieldvalue('toltip',0): 170 175 print "moving tip-domain intersection point" 171 176 … … 179 184 #OK, now we can add our own rift 180 185 nods=rifti['nods']-1 181 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1))))))182 bamg_geometry.Edges=n umpy.vstack((bamg_geometry.Edges,\183 n umpy.array([[pos,count+1,(1+i)]]),\184 n umpy.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))))))186 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1)))))) 187 bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\ 188 np.array([[pos,count+1,(1+i)]]),\ 189 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)))))) 185 190 count+=nods 186 191 … … 189 194 else: 190 195 #Add intersection point to Vertices 191 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,numpy.array([[x,y,1]])))196 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.array([[x,y,1]]))) 192 197 count+=1 193 198 194 199 #Decompose the crossing edge into 2 subedges 195 pos=n umpy.nonzero(numpy.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0]200 pos=np.nonzero(np.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0] 196 201 if not pos: 197 202 raise RuntimeError("bamg error message: a problem occurred...") 198 bamg_geometry.Edges=n umpy.vstack((bamg_geometry.Edges[0:pos-1,:],\199 n umpy.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\200 n umpy.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\203 bamg_geometry.Edges=np.vstack((bamg_geometry.Edges[0:pos-1,:],\ 204 np.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\ 205 np.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\ 201 206 bamg_geometry.Edges[pos+1:,:])) 202 207 203 208 #OK, now we can add our own rift 204 209 nods=rifti['nods']-1 205 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1))))))206 bamg_geometry.Edges=n umpy.vstack((bamg_geometry.Edges,\207 n umpy.array([[count,count+1,2]]),\208 n umpy.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))))))210 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1)))))) 211 bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\ 212 np.array([[count,count+1,2]]),\ 213 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)))))) 209 214 count+=nods 210 215 … … 213 218 else: 214 219 nods=rifti['nods']-1 215 bamg_geometry.Vertices=n umpy.vstack(bamg_geometry.Vertices, numpy.hstack(rifti['x'][:],rifti['y'][:],numpy.ones((nods+1,1))))216 bamg_geometry.Edges =n umpy.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))))220 bamg_geometry.Vertices=np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:],rifti['y'][:],np.ones((nods+1,1)))) 221 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)))) 217 222 count=+nods+1 218 223 … … 224 229 if all(isinstance(track,(str,unicode))): 225 230 A=expread(track) 226 track=n umpy.hstack((A.x.reshape(-1,1),A.y.reshape(-1,1)))231 track=np.hstack((A.x.reshape(-1,),A.y.reshape(-1,))) 227 232 else: 228 233 track=float(track) #for some reason, it is of class "single" 229 if n umpy.size(track,axis=1)==2:230 track=n umpy.hstack((track,3.*numpy.ones((size(track,axis=0),1))))234 if np.size(track,axis=1)==2: 235 track=np.hstack((track,3.*np.ones((size(track,axis=0),1)))) 231 236 232 237 #only keep those inside 233 238 flags=ContourToNodes(track[:,0],track[:,1],domainfile,0)[0] 234 track=track[n umpy.nonzero(flags),:]239 track=track[np.nonzero(flags),:] 235 240 236 241 #Add all points to bamg_geometry 237 nods=n umpy.size(track,axis=0)238 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,track))239 bamg_geometry.Edges =n umpy.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))))))242 nods=np.size(track,axis=0) 243 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,track)) 244 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)))))) 240 245 241 246 #update counter … … 247 252 #recover RequiredVertices 248 253 requiredvertices=options.getfieldvalue('RequiredVertices') #for some reason, it is of class "single" 249 if n umpy.size(requiredvertices,axis=1)==2:250 requiredvertices=n umpy.hstack((requiredvertices,4.*numpy.ones((numpy.size(requiredvertices,axis=0),1))))254 if np.size(requiredvertices,axis=1)==2: 255 requiredvertices=np.hstack((requiredvertices,4.*np.ones((np.size(requiredvertices,axis=0),1)))) 251 256 252 257 #only keep those inside 253 258 flags=ContourToNodes(requiredvertices[:,0],requiredvertices[:,1],domainfile,0)[0] 254 requiredvertices=requiredvertices[n umpy.nonzero(flags)[0],:]259 requiredvertices=requiredvertices[np.nonzero(flags)[0],:] 255 260 #Add all points to bamg_geometry 256 nods=n umpy.size(requiredvertices,axis=0)257 bamg_geometry.Vertices=n umpy.vstack((bamg_geometry.Vertices,requiredvertices))261 nods=np.size(requiredvertices,axis=0) 262 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,requiredvertices)) 258 263 259 264 #update counter … … 275 280 bamg_mesh=bamgmesh(md.private.bamg['mesh'].__dict__) 276 281 else: 277 bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1)))) 278 bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1)))) 282 bamg_mesh.Vertices=np.vstack((md.mesh.x,md.mesh.y,np.ones((md.mesh.numberofvertices)))).T 283 #bamg_mesh.Vertices=np.hstack((md.mesh.x.reshape(-1,),md.mesh.y.reshape(-1,),np.ones((md.mesh.numberofvertices,1)))) 284 bamg_mesh.Triangles=np.hstack((md.mesh.elements,np.ones((md.mesh.numberofelements,1)))) 279 285 280 286 if isinstance(md.rifts.riftstruct,dict): … … 286 292 bamg_options['coeff']=options.getfieldvalue('coeff',1.) 287 293 bamg_options['cutoff']=options.getfieldvalue('cutoff',10.**-5) 288 bamg_options['err']=options.getfieldvalue('err',n umpy.array([[0.01]]))294 bamg_options['err']=options.getfieldvalue('err',np.array([[0.01]])) 289 295 bamg_options['errg']=options.getfieldvalue('errg',0.1) 290 bamg_options['field']=options.getfieldvalue('field',n umpy.empty((0,1)))296 bamg_options['field']=options.getfieldvalue('field',np.empty((0,1))) 291 297 bamg_options['gradation']=options.getfieldvalue('gradation',1.5) 292 298 bamg_options['Hessiantype']=options.getfieldvalue('Hessiantype',0) 293 299 bamg_options['hmin']=options.getfieldvalue('hmin',10.**-100) 294 300 bamg_options['hmax']=options.getfieldvalue('hmax',10.**100) 295 bamg_options['hminVertices']=options.getfieldvalue('hminVertices',n umpy.empty((0,1)))296 bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',n umpy.empty((0,1)))297 bamg_options['hVertices']=options.getfieldvalue('hVertices',n umpy.empty((0,1)))301 bamg_options['hminVertices']=options.getfieldvalue('hminVertices',np.empty((0,1))) 302 bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',np.empty((0,1))) 303 bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1))) 298 304 bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1) 299 305 bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.) 300 306 bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6) 301 307 bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.) 302 bamg_options['metric']=options.getfieldvalue('metric',n umpy.empty((0,1)))308 bamg_options['metric']=options.getfieldvalue('metric',np.empty((0,1))) 303 309 bamg_options['Metrictype']=options.getfieldvalue('Metrictype',0) 304 310 bamg_options['nbjacobi']=options.getfieldvalue('nbjacobi',1) … … 328 334 329 335 #Fill in rest of fields: 330 md.mesh.numberofelements=n umpy.size(md.mesh.elements,axis=0)331 md.mesh.numberofvertices=n umpy.size(md.mesh.x)332 md.mesh.numberofedges=n umpy.size(md.mesh.edges,axis=0)333 md.mesh.vertexonboundary=n umpy.zeros(md.mesh.numberofvertices,bool)336 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 337 md.mesh.numberofvertices=np.size(md.mesh.x) 338 md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) 339 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) 334 340 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 335 341 md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity 336 md.mesh.elementconnectivity[n umpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0342 md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0 337 343 md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) 338 344 339 345 #Check for orphan 340 if n umpy.any(numpy.logical_not(numpy.in1d(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):346 if np.any(np.logical_not(np.in1d(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))): 341 347 raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)") 342 348 … … 349 355 print "Checking Edge crossing..." 350 356 i=0 351 while (i<n umpy.size(geom.Edges,axis=0)):357 while (i<np.size(geom.Edges,axis=0)): 352 358 353 359 #edge counter … … 362 368 363 369 j=i #test edges located AFTER i only 364 while (j<n umpy.size(geom.Edges,axis=0)):370 while (j<np.size(geom.Edges,axis=0)): 365 371 366 372 #edge counter … … 379 385 380 386 #Check if the two edges are crossing one another 381 if SegIntersect(n umpy.array([[x1,y1],[x2,y2]]),numpy.array([[x3,y3],[x4,y4]])):387 if SegIntersect(np.array([[x1,y1],[x2,y2]]),np.array([[x3,y3],[x4,y4]])): 382 388 383 389 #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 384 x=n umpy.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]])))385 y=n umpy.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]])))390 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]]))) 391 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]]))) 386 392 387 393 #Add vertex to the list of vertices 388 geom.Vertices=n umpy.vstack((geom.Vertices,[x,y,min(color1,color2)]))389 id=n umpy.size(geom.Vertices,axis=0)394 geom.Vertices=np.vstack((geom.Vertices,[x,y,min(color1,color2)])) 395 id=np.size(geom.Vertices,axis=0) 390 396 391 397 #Update edges i and j … … 393 399 edgej=geom.Edges[j,:].copy() 394 400 geom.Edges[i,:] =[edgei(0),id ,edgei(2)] 395 geom.Edges=n umpy.vstack((geom.Edges,[id ,edgei(1),edgei(2)]))401 geom.Edges=np.vstack((geom.Edges,[id ,edgei(1),edgei(2)])) 396 402 geom.Edges[j,:] =[edgej(0),id ,edgej(2)] 397 geom.Edges=n umpy.vstack((geom.Edges,[id ,edgej(1),edgej(2)]))403 geom.Edges=np.vstack((geom.Edges,[id ,edgej(1),edgej(2)])) 398 404 399 405 #update current edge second tip … … 405 411 i=0 406 412 num=0 407 while (i<n umpy.size(geom.Vertices,axis=0)):413 while (i<np.size(geom.Vertices,axis=0)): 408 414 409 415 #vertex counter … … 423 429 424 430 #update edges 425 posedges=n umpy.nonzero(geom.Edges==i)431 posedges=np.nonzero(geom.Edges==i) 426 432 geom.Edges[posedges[0],:]=[] 427 posedges=n umpy.nonzero(geom.Edges>i)433 posedges=np.nonzero(geom.Edges>i) 428 434 geom.Edges[posedges]=geom.Edges[posedges]-1 429 435 -
issm/trunk-jpl/src/m/mesh/meshconvert.py
r17558 r21303 1 import numpy 1 import numpy as np 2 2 from collections import OrderedDict 3 3 from BamgConvertMesh import BamgConvertMesh … … 43 43 44 44 #Fill in rest of fields: 45 md.mesh.numberofelements = n umpy.size(md.mesh.elements,axis=0)46 md.mesh.numberofvertices = n umpy.size(md.mesh.x)47 md.mesh.numberofedges = n umpy.size(md.mesh.edges,axis=0)48 md.mesh.vertexonboundary = n umpy.zeros(md.mesh.numberofvertices,bool)45 md.mesh.numberofelements = np.size(md.mesh.elements,axis=0) 46 md.mesh.numberofvertices = np.size(md.mesh.x) 47 md.mesh.numberofedges = np.size(md.mesh.edges,axis=0) 48 md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices,bool) 49 49 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True 50 50 -
issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py
r21069 r21303 1 1 from MatlabFuncs import * 2 2 from model import * 3 from n umpyimport *3 from np.import * 4 4 from pairoptions import * 5 5 from mesh3dsurface import * -
issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py
r17558 r21303 1 import numpy 1 import numpy as np 2 2 from ElementsFromEdge import ElementsFromEdge 3 3 import MatlabFuncs as m … … 19 19 20 20 tips=rift.tips 21 outsidetips=tips[n umpy.nonzero(flags[rift.tips-1])[0]]21 outsidetips=tips[np.nonzero(flags[rift.tips-1])[0]] 22 22 23 23 #we have found outsidetips, tips that touch the domain outline. go through them … … 26 26 #find tip in the segments, take first segment (there should be 2) that holds tip, 27 27 #and node_connected_to_tip is the other node on this segment: 28 tipindex=n umpy.nonzero(rift.segments[:,0]==tip)[0]28 tipindex=np.nonzero(rift.segments[:,0]==tip)[0] 29 29 if tipindex: 30 30 tipindex=tipindex[0] 31 31 node_connected_to_tip=rift.segments[tipindex,1] 32 32 else: 33 tipindex=n umpy.nonzero(rift.segments[:,1]==tip)[0]33 tipindex=np.nonzero(rift.segments[:,1]==tip)[0] 34 34 tipindex=tipindex[0] 35 35 node_connected_to_tip=rift.segments[tipindex,1] … … 41 41 B=node_connected_to_tip 42 42 43 elements=n umpy.empty(0,int)43 elements=np.empty(0,int) 44 44 45 45 while flags(B): #as long as B does not belong to the domain outline, keep looking. … … 48 48 #rule out those we already detected 49 49 already_detected=m.ismember(edgeelements,elements) 50 nextelement=edgeelements(n umpy.nonzero(numpy.logical_not(already_detected))[0])50 nextelement=edgeelements(np.nonzero(np.logical_not(already_detected))[0]) 51 51 #add new detected element to the list of elements we are looking for. 52 elements=n umpy.concatenate((elements,nextelement))52 elements=np.concatenate((elements,nextelement)) 53 53 #new B: 54 B=md.mesh.elements[nextelement-1,n umpy.nonzero(numpy.logical_not(m.ismember(md.mesh.elements[nextelement-1,:],numpy.array([A,B]))))]54 B=md.mesh.elements[nextelement-1,np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement-1,:],np.array([A,B]))))] 55 55 56 56 #take the list of elements on one side of the rift that connect to the tip, 57 57 #and duplicate the tip on them, so as to open the rift to the outside. 58 num=n umpy.size(md.mesh.x)+159 md.mesh.x=n umpy.concatenate((md.mesh.x,md.mesh.x[tip]))60 md.mesh.y=n umpy.concatenate((md.mesh.y,md.mesh.y[tip]))58 num=np.size(md.mesh.x)+1 59 md.mesh.x=np.concatenate((md.mesh.x,md.mesh.x[tip])) 60 md.mesh.y=np.concatenate((md.mesh.y,md.mesh.y[tip])) 61 61 md.mesh.numberofvertices=num 62 62 63 63 #replace tip in elements 64 64 newelements=md.mesh.elements[elements-1,:] 65 pos=n umpy.nonzero(newelements==tip)65 pos=np.nonzero(newelements==tip) 66 66 newelements[pos]=num 67 67 md.mesh.elements[elements-1,:]=newelements 68 rift.tips=n umpy.concatenate((rift.tips,num))68 rift.tips=np.concatenate((rift.tips,num)) 69 69 70 70 #deal with segments 71 tipsegments=n umpy.nonzero(numpy.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0]71 tipsegments=np.nonzero(np.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0] 72 72 for segment_index in tipsegments: 73 pos=n umpy.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0]73 pos=np.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0] 74 74 other_node=md.mesh.segments[segment_index,pos] 75 75 if not isconnected(md.mesh.elements,other_node,tip): 76 pos=n umpy.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0]76 pos=np.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0] 77 77 md.mesh.segments[segment_index,pos]=num 78 78 79 79 #Fill in rest of fields: 80 md.mesh.numberofelements=n umpy.size(md.mesh.elements,axis=0)81 md.mesh.numberofvertices=n umpy.size(md.mesh.x)82 md.mesh.vertexonboundary=n umpy.zeros(numpy.size(md.mesh.x),bool)80 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 81 md.mesh.numberofvertices=np.size(md.mesh.x) 82 md.mesh.vertexonboundary=np.zeros(np.size(md.mesh.x),bool) 83 83 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 84 84 md.rifts.numrifts=length(md.rifts.riftstruct) -
issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py
r20910 r21303 1 import numpy 1 import numpy as np 2 2 from TriMeshProcessRifts import TriMeshProcessRifts 3 3 from ContourToMesh import ContourToMesh … … 33 33 #Fill in rest of fields: 34 34 numrifts=len(md.rifts.riftstruct) 35 md.mesh.numberofelements=n umpy.size(md.mesh.elements,axis=0)36 md.mesh.numberofvertices=n umpy.size(md.mesh.x)37 md.mesh.vertexonboundary=n umpy.zeros(numpy.size(md.mesh.x),bool)35 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 36 md.mesh.numberofvertices=np.size(md.mesh.x) 37 md.mesh.vertexonboundary=np.zeros(np.size(md.mesh.x),bool) 38 38 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 39 39 40 40 #get coordinates of rift tips 41 41 for rift in md.rifts.riftstruct: 42 rift['tip1coordinates']=n umpy.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)))43 rift['tip2coordinates']=n umpy.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)))42 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,))) 43 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,))) 44 44 45 45 #In case we have rifts that open up the domain outline, we need to open them: … … 58 58 #get elements that are not correctly oriented in the correct direction: 59 59 aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y) 60 pos=n umpy.nonzero(aires<0)[0]61 md.mesh.elements[pos,:]=n umpy.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)))60 pos=np.nonzero(aires<0)[0] 61 md.mesh.elements[pos,:]=np.vstack((md.mesh.elements[pos,1],md.mesh.elements[pos,0],md.mesh.elements[pos,2])).T 62 62 63 63 return md -
issm/trunk-jpl/src/m/mesh/roundmesh.py
r19433 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 from collections import OrderedDict … … 20 20 21 21 #Get number of points on the circle 22 pointsonedge=n umpy.floor((2.*numpy.pi*radius) / resolution)22 pointsonedge=np.floor((2.*np.pi*radius) / resolution) 23 23 24 24 #Calculate the cartesians coordinates of the points 25 x_list=n umpy.ones(pointsonedge)26 y_list=n umpy.ones(pointsonedge)27 theta=n umpy.linspace(0.,2.*numpy.pi,num=pointsonedge,endpoint=False)28 x_list=roundsigfig(radius*x_list*n umpy.cos(theta),12)29 y_list=roundsigfig(radius*y_list*n umpy.sin(theta),12)25 x_list=np.ones(pointsonedge) 26 y_list=np.ones(pointsonedge) 27 theta=np.linspace(0.,2.*np.pi,num=pointsonedge,endpoint=False) 28 x_list=roundsigfig(radius*x_list*np.cos(theta),12) 29 y_list=roundsigfig(radius*y_list*np.sin(theta),12) 30 30 A=OrderedDict() 31 31 A['x']=[x_list] … … 39 39 40 40 #move the closest node to the center 41 pos=n umpy.argmin(md.mesh.x**2+md.mesh.y**2)41 pos=np.argmin(md.mesh.x**2+md.mesh.y**2) 42 42 md.mesh.x[pos]=0. 43 43 md.mesh.y[pos]=0. … … 50 50 def roundsigfig(x,n): 51 51 52 digits=n umpy.ceil(numpy.log10(numpy.abs(x)))52 digits=np.ceil(np.log10(np.abs(x))) 53 53 x=x/10.**digits 54 x=n umpy.round(x,decimals=n)54 x=np.round(x,decimals=n) 55 55 x=x*10.**digits 56 56 57 pos=n umpy.nonzero(numpy.isnan(x))57 pos=np.nonzero(np.isnan(x)) 58 58 x[pos]=0. 59 59 -
issm/trunk-jpl/src/m/mesh/squaremesh.py
r20910 r21303 1 import numpy 1 import numpy as np 2 2 from NodeConnectivity import NodeConnectivity 3 3 from ElementConnectivity import ElementConnectivity … … 22 22 23 23 #initialization 24 index=n umpy.zeros((nel,3),int)25 x=n umpy.zeros((nx*ny))26 y=n umpy.zeros((nx*ny))24 index=np.zeros((nel,3),int) 25 x=np.zeros((nx*ny)) 26 y=np.zeros((nx*ny)) 27 27 28 28 #create coordinates … … 43 43 44 44 #Scale x and y 45 x=x/n umpy.max(x)*Lx46 y=y/n umpy.max(y)*Ly45 x=x/np.max(x)*Lx 46 y=y/np.max(y)*Ly 47 47 48 48 #create segments 49 segments=n umpy.zeros((2*(nx-1)+2*(ny-1),3),int)49 segments=np.zeros((2*(nx-1)+2*(ny-1),3),int) 50 50 #left edge: 51 segments[0:ny-1,:]=n umpy.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)))51 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,))) 52 52 #right edge: 53 segments[ny-1:2*(ny-1),:]=n umpy.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)))53 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,))) 54 54 #front edge: 55 segments[2*(ny-1):2*(ny-1)+(nx-1),:]=n umpy.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)))55 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,))) 56 56 #back edge 57 segments[2*(ny-1)+(nx-1):2*(nx-1)+2*(ny-1),:]=n umpy.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)))57 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,))) 58 58 59 59 #plug coordinates and nodes … … 62 62 md.mesh.y=y 63 63 md.mesh.numberofvertices=nods 64 md.mesh.vertexonboundary=n umpy.zeros((nods),bool)64 md.mesh.vertexonboundary=np.zeros((nods),bool) 65 65 md.mesh.vertexonboundary[segments[:,0:2]-1]=True 66 66 -
issm/trunk-jpl/src/m/mesh/triangle.py
r21174 r21303 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 from mesh2d import mesh2d 4 4 from TriMesh import TriMesh … … 56 56 57 57 #Fill in rest of fields: 58 md.mesh.numberofelements = n umpy.size(md.mesh.elements,axis=0)59 md.mesh.numberofvertices = n umpy.size(md.mesh.x)60 md.mesh.vertexonboundary = n umpy.zeros(md.mesh.numberofvertices,bool)58 md.mesh.numberofelements = np.size(md.mesh.elements,axis=0) 59 md.mesh.numberofvertices = np.size(md.mesh.x) 60 md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices,bool) 61 61 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True 62 62 -
issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py
r17264 r21303 49 49 50 50 def ismember(a,s): 51 import numpy 51 import numpy as np 52 52 53 if not isinstance(s,(tuple,list,dict,n umpy.ndarray)):53 if not isinstance(s,(tuple,list,dict,np.ndarray)): 54 54 s=[s] 55 55 56 if not isinstance(a,(tuple,list,dict,n umpy.ndarray)):56 if not isinstance(a,(tuple,list,dict,np.ndarray)): 57 57 a=[a] 58 58 59 if not isinstance(a,n umpy.ndarray):59 if not isinstance(a,np.ndarray): 60 60 b=[item in s for item in a] 61 61 62 62 else: 63 if not isinstance(s,n umpy.ndarray):64 b=n umpy.empty_like(a)63 if not isinstance(s,np.ndarray): 64 b=np.empty_like(a) 65 65 for i,item in enumerate(a.flat): 66 66 b.flat[i]=item in s 67 67 else: 68 b=n umpy.in1d(a.flat,s.flat).reshape(a.shape)68 b=np.in1d(a.flat,s.flat).reshape(a.shape) 69 69 70 70 return b 71 71 72 72 def det(a): 73 import numpy 73 import numpy as np 74 74 75 75 if a.shape==(1,): … … 83 83 84 84 def sparse(ivec,jvec,svec,m=0,n=0,nzmax=0): 85 import numpy 85 import numpy as np 86 86 87 87 if not m: 88 m=n umpy.max(ivec)88 m=np.max(ivec) 89 89 if not n: 90 n=n umpy.max(jvec)90 n=np.max(jvec) 91 91 92 a=n umpy.zeros((m,n))92 a=np.zeros((m,n)) 93 93 94 94 for i,j,s in zip(ivec.reshape(-1,order='F'),jvec.reshape(-1,order='F'),svec.reshape(-1,order='F')): … … 98 98 99 99 def heaviside(x): 100 import numpy 100 import numpy as np 101 101 102 y=n umpy.zeros_like(x)103 y[n umpy.nonzero(x> 0.)]=1.104 y[n umpy.nonzero(x==0.)]=0.5102 y=np.zeros_like(x) 103 y[np.nonzero(x> 0.)]=1. 104 y[np.nonzero(x==0.)]=0.5 105 105 106 106 return y -
issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py
r14211 r21303 1 import numpy as np 2 1 3 def logical_and_n(*arg): 2 from numpy import logical_and 3 4 4 5 if len(arg): 5 6 result=arg[0] … … 12 13 13 14 def logical_or_n(*arg): 14 from numpy import logical_or 15 15 16 16 if len(arg): 17 17 result=arg[0] -
issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py
r17480 r21303 1 1 #Module import 2 import numpy 2 import numpy as np 3 3 from math import isnan 4 4 import MatlabFuncs as m … … 33 33 34 34 #matrix 35 elif isinstance(field,n umpy.ndarray):35 elif isinstance(field,np.ndarray): 36 36 string=displayunit(offset,name,str(field.shape),comment) 37 37 -
issm/trunk-jpl/src/m/miscellaneous/isnans.py
r13011 r21303 1 import numpy 1 import numpy as np 2 2 3 3 def isnans(array): … … 13 13 returnvalue=0 14 14 else: 15 returnvalue=n umpy.isnan(array)15 returnvalue=np.isnan(array) 16 16 17 17 return returnvalue -
issm/trunk-jpl/src/m/parameterization/contourenvelope.py
r20910 r21303 1 1 import os.path 2 import numpy 2 import numpy as np 3 3 import copy 4 4 from NodeConnectivity import NodeConnectivity … … 40 40 #Now, build the connectivity tables for this mesh. 41 41 #Computing connectivity 42 if n umpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d:42 if np.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d: 43 43 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices)[0] 44 if n umpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d:44 if np.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and np.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d: 45 45 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity)[0] 46 46 … … 65 65 #get flag list of elements and nodes inside the contour 66 66 nodein=ContourToMesh(elements,x,y,file,'node',1) 67 elemin=(n umpy.sum(nodein(elements),axis=1)==numpy.size(elements,axis=1))67 elemin=(np.sum(nodein(elements),axis=1)==np.size(elements,axis=1)) 68 68 #modify element connectivity 69 elemout=n umpy.nonzero(numpy.logical_not(elemin))[0]69 elemout=np.nonzero(np.logical_not(elemin))[0] 70 70 elementconnectivity[elemout,:]=0 71 elementconnectivity[n umpy.nonzero(m.ismember(elementconnectivity,elemout+1))]=071 elementconnectivity[np.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 72 72 else: 73 73 #get flag list of elements and nodes inside the contour 74 nodein=n umpy.zeros(numberofvertices)75 elemin=n umpy.zeros(numberofelements)74 nodein=np.zeros(numberofvertices) 75 elemin=np.zeros(numberofelements) 76 76 77 pos=n umpy.nonzero(flags)77 pos=np.nonzero(flags) 78 78 elemin[pos]=1 79 79 nodein[elements[pos,:]-1]=1 80 80 81 81 #modify element connectivity 82 elemout=n umpy.nonzero(numpy.logical_not(elemin))[0]82 elemout=np.nonzero(np.logical_not(elemin))[0] 83 83 elementconnectivity[elemout,:]=0 84 elementconnectivity[n umpy.nonzero(m.ismember(elementconnectivity,elemout+1))]=084 elementconnectivity[np.nonzero(m.ismember(elementconnectivity,elemout+1))]=0 85 85 86 86 #Find element on boundary … … 88 88 flag=copy.deepcopy(elementconnectivity) 89 89 if len(args)==1: 90 flag[n umpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]]91 elementonboundary=n umpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0)90 flag[np.nonzero(flag)]=elemin[flag[np.nonzero(flag)]] 91 elementonboundary=np.logical_and(np.prod(flag,axis=1)==0,np.sum(flag,axis=1)>0) 92 92 93 93 #Find segments on boundary 94 pos=n umpy.nonzero(elementonboundary)[0]95 num_segments=n umpy.size(pos)96 segments=n umpy.zeros((num_segments*3,3),int)94 pos=np.nonzero(elementonboundary)[0] 95 num_segments=np.size(pos) 96 segments=np.zeros((num_segments*3,3),int) 97 97 count=0 98 98 99 99 for el1 in pos: 100 els2=elementconnectivity[el1,n umpy.nonzero(elementconnectivity[el1,:])[0]]-1101 if n umpy.size(els2)>1:102 flag=n umpy.intersect1d(numpy.intersect1d(elements[els2[0],:],elements[els2[1],:]),elements[el1,:])100 els2=elementconnectivity[el1,np.nonzero(elementconnectivity[el1,:])[0]]-1 101 if np.size(els2)>1: 102 flag=np.intersect1d(np.intersect1d(elements[els2[0],:],elements[els2[1],:]),elements[el1,:]) 103 103 nods1=elements[el1,:] 104 nods1=n umpy.delete(nods1,numpy.nonzero(nods1==flag))104 nods1=np.delete(nods1,np.nonzero(nods1==flag)) 105 105 segments[count,:]=[nods1[0],nods1[1],el1+1] 106 106 107 ord1=n umpy.nonzero(nods1[0]==elements[el1,:])[0][0]108 ord2=n umpy.nonzero(nods1[1]==elements[el1,:])[0][0]107 ord1=np.nonzero(nods1[0]==elements[el1,:])[0][0] 108 ord2=np.nonzero(nods1[1]==elements[el1,:])[0][0] 109 109 110 110 #swap segment nodes if necessary … … 113 113 segments[count,0]=segments[count,1] 114 114 segments[count,1]=temp 115 segments[count,0:2]=n umpy.flipud(segments[count,0:2])115 segments[count,0:2]=np.flipud(segments[count,0:2]) 116 116 count+=1 117 117 else: 118 118 nods1=elements[el1,:] 119 flag=n umpy.setdiff1d(nods1,elements[els2,:])119 flag=np.setdiff1d(nods1,elements[els2,:]) 120 120 for j in xrange(0,3): 121 nods=n umpy.delete(nods1,j)122 if n umpy.any(m.ismember(flag,nods)):121 nods=np.delete(nods1,j) 122 if np.any(m.ismember(flag,nods)): 123 123 segments[count,:]=[nods[0],nods[1],el1+1] 124 ord1=n umpy.nonzero(nods[0]==elements[el1,:])[0][0]125 ord2=n umpy.nonzero(nods[1]==elements[el1,:])[0][0]124 ord1=np.nonzero(nods[0]==elements[el1,:])[0][0] 125 ord2=np.nonzero(nods[1]==elements[el1,:])[0][0] 126 126 if ( (ord1==0 and ord2==1) or (ord1==1 and ord2==2) or (ord1==2 and ord2==0) ): 127 127 temp=segments[count,0] 128 128 segments[count,0]=segments[count,1] 129 129 segments[count,1]=temp 130 segments[count,0:2]=n umpy.flipud(segments[count,0:2])130 segments[count,0:2]=np.flipud(segments[count,0:2]) 131 131 count+=1 132 132 segments=segments[0:count,:] -
issm/trunk-jpl/src/m/parameterization/setflowequation.py
r19535 r21303 1 import numpy 1 import numpy as np 2 2 from model import model 3 3 from pairoptions import pairoptions 4 import MatlabFuncs as m5 import PythonFuncs as p6 4 from FlagElements import FlagElements 7 5 … … 33 31 #process options 34 32 options=pairoptions(*args) 35 # options=deleteduplicates(options,1);33 # options=deleteduplicates(options,1); 36 34 37 35 #Find_out what kind of coupling to use 38 36 coupling_method=options.getfieldvalue('coupling','tiling') 39 if not m.strcmpi(coupling_method,'tiling') and not m.strcmpi(coupling_method,'penalties'):37 if not coupling_method in ['tiling','penalties']: 40 38 raise TypeError("coupling type can only be: tiling or penalties") 41 39 … … 47 45 FSflag = FlagElements(md,options.getfieldvalue('FS','')) 48 46 filltype = options.getfieldvalue('fill','none') 49 50 47 #Flag the elements that have not been flagged as filltype 51 if m.strcmpi(filltype,'SIA'): 52 SIAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SSAflag,HOflag)))]=True 53 elif m.strcmpi(filltype,'SSA'): 54 SSAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,HOflag,FSflag)))]=True 55 elif m.strcmpi(filltype,'HO'): 56 HOflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,SSAflag,FSflag)))]=True 57 48 if 'SIA' in filltype: 49 SIAflag= ~SSAflag & ~HOflag 50 elif 'SSA' in filltype: 51 SSAflag=~SIAflag & ~HOflag & ~FSflag 52 elif 'HO' in filltype: 53 HOflag=~SIAflag & ~SSAflag & ~FSflag 58 54 #check that each element has at least one flag 59 55 if not any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag): … … 63 59 if any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag>1): 64 60 print "setflowequation warning message: some elements have several types, higher order type is used for them" 65 SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,SSAflag))]=False 66 SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,HOflag))]=False 67 SSAflag[numpy.nonzero(numpy.logical_and(SSAflag,HOflag))]=False 61 SIAflag[np.where(np.logical_and(SIAflag,SSAflag))]=False 62 SIAflag[np.where(np.logical_and(SIAflag,HOflag))]=False 63 SSAflag[np.where(np.logical_and(SSAflag,HOflag))]=False 64 65 #check that L1L2 is not coupled to any other model for now 66 if any(L1L2flag) and any(SIAflag+SSAflag+HOflag+FSflag): 67 raise TypeError('L1L2 cannot be coupled to any other model') 68 69 #Check that no HO or FS for 2d mesh 70 if domaintype(md.mesh)=='2Dhorizontal': 71 if any(FSflag+HOflag): 72 raise TypeError('FS and HO elements not allowed in 2d mesh, extrude it first') 68 73 69 74 #FS can only be used alone for now: … … 72 77 73 78 #Initialize node fields 74 nodeonSIA=n umpy.zeros(md.mesh.numberofvertices,bool)75 nodeonSIA[md.mesh.elements[n umpy.nonzero(SIAflag),:]-1]=True76 nodeonSSA=n umpy.zeros(md.mesh.numberofvertices,bool)77 nodeonSSA[md.mesh.elements[n umpy.nonzero(SSAflag),:]-1]=True78 nodeonL1L2=n umpy.zeros(md.mesh.numberofvertices,bool)79 nodeonL1L2[md.mesh.elements[n umpy.nonzero(L1L2flag),:]-1]=True80 nodeonHO=n umpy.zeros(md.mesh.numberofvertices,bool)81 nodeonHO[md.mesh.elements[n umpy.nonzero(HOflag),:]-1]=True82 nodeonFS=n umpy.zeros(md.mesh.numberofvertices,bool)83 noneflag=n umpy.zeros(md.mesh.numberofelements,bool)79 nodeonSIA=np.zeros(md.mesh.numberofvertices,bool) 80 nodeonSIA[md.mesh.elements[np.where(SIAflag),:]-1]=True 81 nodeonSSA=np.zeros(md.mesh.numberofvertices,bool) 82 nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True 83 nodeonL1L2=np.zeros(md.mesh.numberofvertices,bool) 84 nodeonL1L2[md.mesh.elements[np.where(L1L2flag),:]-1]=True 85 nodeonHO=np.zeros(md.mesh.numberofvertices,bool) 86 nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True 87 nodeonFS=np.zeros(md.mesh.numberofvertices,bool) 88 noneflag=np.zeros(md.mesh.numberofelements,bool) 84 89 85 90 #First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA) 86 91 if any(FSflag): 87 # 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 88 fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)).astype(int)+ \ 89 numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)).astype(int)+ \ 90 numpy.logical_not(numpy.isnan(md.stressbalance.spcvz)).astype(int)==3, \ 91 numpy.logical_and(nodeonHO,nodeonFS)).astype(int) #find all the nodes on the boundary of the domain without icefront 92 # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront 93 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 94 FSflag[numpy.nonzero(fullspcelems.reshape(-1))]=False 95 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True 92 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 93 fullspcelems=np.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6 #find all the nodes on the boundary of the domain without icefront 94 FSflag[np.where(fullspcelems.reshape(-1))]=False 95 nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True 96 96 97 97 #Then complete with NoneApproximation or the other model used if there is no FS 98 98 if any(FSflag): 99 99 if any(HOflag): #fill with HO 100 HOflag[ numpy.logical_not(FSflag)]=True101 nodeonHO[md.mesh.elements[n umpy.nonzero(HOflag),:]-1]=True100 HOflag[~FSflag]=True 101 nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True 102 102 elif any(SSAflag): #fill with SSA 103 SSAflag[ numpy.logical_not(FSflag)]=True104 nodeonSSA[md.mesh.elements[n umpy.nonzero(SSAflag),:]-1]=True103 SSAflag[~FSflag]=True 104 nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True 105 105 else: #fill with none 106 noneflag[n umpy.nonzero(numpy.logical_not(FSflag))]=True106 noneflag[np.where(~FSflag)]=True 107 107 108 108 #Now take care of the coupling between SSA and HO 109 md.stressbalance.vertex_pairing=n umpy.array([])110 nodeonSSAHO=n umpy.zeros(md.mesh.numberofvertices,bool)111 nodeonHOFS=n umpy.zeros(md.mesh.numberofvertices,bool)112 nodeonSSAFS=n umpy.zeros(md.mesh.numberofvertices,bool)113 SSAHOflag=n umpy.zeros(md.mesh.numberofelements,bool)114 SSAFSflag=n umpy.zeros(md.mesh.numberofelements,bool)115 HOFSflag=n umpy.zeros(md.mesh.numberofelements,bool)116 if m.strcmpi(coupling_method,'penalties'):109 md.stressbalance.vertex_pairing=np.array([]) 110 nodeonSSAHO=np.zeros(md.mesh.numberofvertices,bool) 111 nodeonHOFS=np.zeros(md.mesh.numberofvertices,bool) 112 nodeonSSAFS=np.zeros(md.mesh.numberofvertices,bool) 113 SSAHOflag=np.zeros(md.mesh.numberofelements,bool) 114 SSAFSflag=np.zeros(md.mesh.numberofelements,bool) 115 HOFSflag=np.zeros(md.mesh.numberofelements,bool) 116 if coupling_method=='penalties': 117 117 #Create the border nodes between HO and SSA and extrude them 118 118 numnodes2d=md.mesh.numberofvertices2d 119 119 numlayers=md.mesh.numberoflayers 120 bordernodes2d=n umpy.nonzero(numpy.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1 #Nodes connected to two different types of elements120 bordernodes2d=np.where(np.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1 #Nodes connected to two different types of elements 121 121 122 122 #initialize and fill in penalties structure 123 if n umpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))):124 penalties=n umpy.zeros((0,2))123 if np.all(np.logical_not(np.isnan(bordernodes2d))): 124 penalties=np.zeros((0,2)) 125 125 for i in xrange(1,numlayers): 126 penalties=n umpy.vstack((penalties,numpy.hstack((bordernodes2d.reshape(-1,1),bordernodes2d.reshape(-1,1)+md.mesh.numberofvertices2d*(i)))))126 penalties=np.vstack((penalties,np.vstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i))).T)) 127 127 md.stressbalance.vertex_pairing=penalties 128 128 129 elif m.strcmpi(coupling_method,'tiling'):130 if 129 elif coupling_method=='tiling': 130 if any(SSAflag) and any(HOflag): #coupling SSA HO 131 131 #Find node at the border 132 nodeonSSAHO[n umpy.nonzero(numpy.logical_and(nodeonSSA,nodeonHO))]=True132 nodeonSSAHO[np.where(np.logical_and(nodeonSSA,nodeonHO))]=True 133 133 #SSA elements in contact with this layer become SSAHO elements 134 matrixelements= m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAHO)[0])135 commonelements=n umpy.sum(matrixelements,axis=1)!=0136 commonelements[n umpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA137 SSAflag[n umpy.nonzero(commonelements)]=False #these elements are now SSAHOelements138 SSAHOflag[n umpy.nonzero(commonelements)]=True134 matrixelements=nodeonSSAHO[md.mesh.elements-1] 135 commonelements=np.sum(matrixelements,axis=1)!=0 136 commonelements[np.where(HOflag)]=False #only one layer: the elements previously in SSA 137 SSAflag[np.where(commonelements)]=False #these elements are now SSAHOelements 138 SSAHOflag[np.where(commonelements)]=True 139 139 nodeonSSA[:]=False 140 nodeonSSA[md.mesh.elements[n umpy.nonzero(SSAflag),:]-1]=True140 nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True 141 141 142 142 #rule out elements that don't touch the 2 boundaries 143 pos=n umpy.nonzero(SSAHOflag)[0]144 elist=n umpy.zeros(numpy.size(pos),dtype=int)145 elist = elist + n umpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)146 elist = elist - n umpy.sum(nodeonHO[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool)147 pos1=n umpy.nonzero(elist==1)[0]143 pos=np.where(SSAHOflag)[0] 144 elist=np.zeros(np.size(pos),dtype=int) 145 elist = elist + np.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 146 elist = elist - np.sum(nodeonHO[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) 147 pos1=np.where(elist==1)[0] 148 148 SSAflag[pos[pos1]]=True 149 149 SSAHOflag[pos[pos1]]=False 150 pos2=n umpy.nonzero(elist==-1)[0]150 pos2=np.where(elist==-1)[0] 151 151 HOflag[pos[pos2]]=True 152 152 SSAHOflag[pos[pos2]]=False … … 154 154 #Recompute nodes associated to these elements 155 155 nodeonSSA[:]=False 156 nodeonSSA[md.mesh.elements[n umpy.nonzero(SSAflag),:]-1]=True156 nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True 157 157 nodeonHO[:]=False 158 nodeonHO[md.mesh.elements[n umpy.nonzero(HOflag),:]-1]=True158 nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True 159 159 nodeonSSAHO[:]=False 160 nodeonSSAHO[md.mesh.elements[n umpy.nonzero(SSAHOflag),:]-1]=True160 nodeonSSAHO[md.mesh.elements[np.where(SSAHOflag),:]-1]=True 161 161 162 162 elif any(HOflag) and any(FSflag): #coupling HO FS 163 163 #Find node at the border 164 nodeonHOFS[n umpy.nonzero(numpy.logical_and(nodeonHO,nodeonFS))]=True164 nodeonHOFS[np.where(np.logical_and(nodeonHO,nodeonFS))]=True 165 165 #FS elements in contact with this layer become HOFS elements 166 matrixelements= m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonHOFS)[0])167 commonelements=n umpy.sum(matrixelements,axis=1)!=0168 commonelements[n umpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA169 FSflag[n umpy.nonzero(commonelements)]=False #these elements are now SSAHOelements170 HOFSflag[n umpy.nonzero(commonelements)]=True171 nodeonFS=n umpy.zeros(md.mesh.numberofvertices,bool)172 nodeonFS[md.mesh.elements[n umpy.nonzero(FSflag),:]-1]=True166 matrixelements=nodeonHOFS[md.mesh.elements-1] 167 commonelements=np.sum(matrixelements,axis=1)!=0 168 commonelements[np.where(HOflag)]=False #only one layer: the elements previously in SSA 169 FSflag[np.where(commonelements)]=False #these elements are now SSAHOelements 170 HOFSflag[np.where(commonelements)]=True 171 nodeonFS=np.zeros(md.mesh.numberofvertices,bool) 172 nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True 173 173 174 174 #rule out elements that don't touch the 2 boundaries 175 pos=n umpy.nonzero(HOFSflag)[0]176 elist=n umpy.zeros(numpy.size(pos),dtype=int)177 elist = elist + n umpy.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool)178 elist = elist - n umpy.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool)179 pos1=n umpy.nonzero(elist==1)[0]175 pos=np.where(HOFSflag)[0] 176 elist=np.zeros(np.size(pos),dtype=int) 177 elist = elist + np.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 178 elist = elist - np.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 179 pos1=np.where(elist==1)[0] 180 180 FSflag[pos[pos1]]=True 181 181 HOFSflag[pos[pos1]]=False 182 pos2=n umpy.nonzero(elist==-1)[0]182 pos2=np.where(elist==-1)[0] 183 183 HOflag[pos[pos2]]=True 184 184 HOFSflag[pos[pos2]]=False … … 186 186 #Recompute nodes associated to these elements 187 187 nodeonFS[:]=False 188 nodeonFS[md.mesh.elements[n umpy.nonzero(FSflag),:]-1]=True188 nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True 189 189 nodeonHO[:]=False 190 nodeonHO[md.mesh.elements[n umpy.nonzero(HOflag),:]-1]=True190 nodeonHO[md.mesh.elements[np.where(HOflag),:]-1]=True 191 191 nodeonHOFS[:]=False 192 nodeonHOFS[md.mesh.elements[numpy.nonzero(HOFSflag),:]-1]=True 193 192 nodeonHOFS[md.mesh.elements[np.where(HOFSflag),:]-1]=True 194 193 elif any(FSflag) and any(SSAflag): 195 194 #Find node at the border 196 nodeonSSAFS[n umpy.nonzero(numpy.logical_and(nodeonSSA,nodeonFS))]=True195 nodeonSSAFS[np.where(np.logical_and(nodeonSSA,nodeonFS))]=True 197 196 #FS elements in contact with this layer become SSAFS elements 198 matrixelements= m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAFS)[0])199 commonelements=n umpy.sum(matrixelements,axis=1)!=0200 commonelements[n umpy.nonzero(SSAflag)]=False #only one layer: the elements previously in SSA201 FSflag[n umpy.nonzero(commonelements)]=False #these elements are now SSASSAelements202 SSAFSflag[n umpy.nonzero(commonelements)]=True203 nodeonFS=n umpy.zeros(md.mesh.numberofvertices,bool)204 nodeonFS[md.mesh.elements[n umpy.nonzero(FSflag),:]-1]=True197 matrixelements=nodeonSSAFS[md.mesh.elements-1] 198 commonelements=np.sum(matrixelements,axis=1)!=0 199 commonelements[np.where(SSAflag)]=False #only one layer: the elements previously in SSA 200 FSflag[np.where(commonelements)]=False #these elements are now SSASSAelements 201 SSAFSflag[np.where(commonelements)]=True 202 nodeonFS=np.zeros(md.mesh.numberofvertices,bool) 203 nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True 205 204 206 205 #rule out elements that don't touch the 2 boundaries 207 pos=n umpy.nonzero(SSAFSflag)[0]208 elist=n umpy.zeros(numpy.size(pos),dtype=int)209 elist = elist + n umpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)210 elist = elist - n umpy.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool)211 pos1=n umpy.nonzero(elist==1)[0]206 pos=np.where(SSAFSflag)[0] 207 elist=np.zeros(np.size(pos),dtype=int) 208 elist = elist + np.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 209 elist = elist - np.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 210 pos1=np.where(elist==1)[0] 212 211 SSAflag[pos[pos1]]=True 213 212 SSAFSflag[pos[pos1]]=False 214 pos2=n umpy.nonzero(elist==-1)[0]213 pos2=np.where(elist==-1)[0] 215 214 FSflag[pos[pos2]]=True 216 215 SSAFSflag[pos[pos2]]=False … … 218 217 #Recompute nodes associated to these elements 219 218 nodeonSSA[:]=False 220 nodeonSSA[md.mesh.elements[n umpy.nonzero(SSAflag),:]-1]=True219 nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True 221 220 nodeonFS[:]=False 222 nodeonFS[md.mesh.elements[n umpy.nonzero(FSflag),:]-1]=True221 nodeonFS[md.mesh.elements[np.where(FSflag),:]-1]=True 223 222 nodeonSSAFS[:]=False 224 nodeonSSAFS[md.mesh.elements[n umpy.nonzero(SSAFSflag),:]-1]=True223 nodeonSSAFS[md.mesh.elements[np.where(SSAFSflag),:]-1]=True 225 224 226 225 elif any(FSflag) and any(SIAflag): … … 228 227 229 228 #Create SSAHOApproximation where needed 230 md.flowequation.element_equation=n umpy.zeros(md.mesh.numberofelements,int)231 md.flowequation.element_equation[n umpy.nonzero(noneflag)]=0232 md.flowequation.element_equation[n umpy.nonzero(SIAflag)]=1233 md.flowequation.element_equation[n umpy.nonzero(SSAflag)]=2234 md.flowequation.element_equation[n umpy.nonzero(L1L2flag)]=3235 md.flowequation.element_equation[n umpy.nonzero(HOflag)]=4236 md.flowequation.element_equation[n umpy.nonzero(FSflag)]=5237 md.flowequation.element_equation[n umpy.nonzero(SSAHOflag)]=6238 md.flowequation.element_equation[n umpy.nonzero(SSAFSflag)]=7239 md.flowequation.element_equation[n umpy.nonzero(HOFSflag)]=8229 md.flowequation.element_equation=np.zeros(md.mesh.numberofelements,int) 230 md.flowequation.element_equation[np.where(noneflag)]=0 231 md.flowequation.element_equation[np.where(SIAflag)]=1 232 md.flowequation.element_equation[np.where(SSAflag)]=2 233 md.flowequation.element_equation[np.where(L1L2flag)]=3 234 md.flowequation.element_equation[np.where(HOflag)]=4 235 md.flowequation.element_equation[np.where(FSflag)]=5 236 md.flowequation.element_equation[np.where(SSAHOflag)]=6 237 md.flowequation.element_equation[np.where(SSAFSflag)]=7 238 md.flowequation.element_equation[np.where(HOFSflag)]=8 240 239 241 240 #border … … 245 244 246 245 #Create vertices_type 247 md.flowequation.vertex_equation=n umpy.zeros(md.mesh.numberofvertices,int)248 pos=n umpy.nonzero(nodeonSSA)246 md.flowequation.vertex_equation=np.zeros(md.mesh.numberofvertices,int) 247 pos=np.where(nodeonSSA) 249 248 md.flowequation.vertex_equation[pos]=2 250 pos=n umpy.nonzero(nodeonL1L2)249 pos=np.where(nodeonL1L2) 251 250 md.flowequation.vertex_equation[pos]=3 252 pos=n umpy.nonzero(nodeonHO)251 pos=np.where(nodeonHO) 253 252 md.flowequation.vertex_equation[pos]=4 254 pos=n umpy.nonzero(nodeonFS)253 pos=np.where(nodeonFS) 255 254 md.flowequation.vertex_equation[pos]=5 256 255 #DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority) 257 pos=n umpy.nonzero(nodeonSIA)256 pos=np.where(nodeonSIA) 258 257 md.flowequation.vertex_equation[pos]=1 259 258 if any(FSflag): 260 pos=n umpy.nonzero(numpy.logical_not(nodeonFS))259 pos=np.where(np.logical_not(nodeonFS)) 261 260 if not (any(HOflag) or any(SSAflag)): 262 261 md.flowequation.vertex_equation[pos]=0 263 pos=n umpy.nonzero(nodeonSSAHO)262 pos=np.where(nodeonSSAHO) 264 263 md.flowequation.vertex_equation[pos]=6 265 pos=n umpy.nonzero(nodeonHOFS)264 pos=np.where(nodeonHOFS) 266 265 md.flowequation.vertex_equation[pos]=7 267 pos=n umpy.nonzero(nodeonSSAFS)266 pos=np.where(nodeonSSAFS) 268 267 md.flowequation.vertex_equation[pos]=8 269 268 -
issm/trunk-jpl/src/m/parameterization/sethydrostaticmask.py
r18294 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 from model import model … … 20 20 """ 21 21 22 if n umpy.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:22 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: 23 23 raise IOError("hydrostaticmask error message: fields in md.geometry do not have the right size.") 24 24 … … 27 27 28 28 #Check consistency of geometry 29 if any(md.geometry.base[n umpy.nonzero(md.mask.groundedice_levelset>0.)]!=md.geometry.bed[numpy.nonzero(md.mask.groundedice_levelset>0.)]):29 if any(md.geometry.base[np.nonzero(md.mask.groundedice_levelset>0.)]!=md.geometry.bed[np.nonzero(md.mask.groundedice_levelset>0.)]): 30 30 print "WARNING: md.geometry.bed and md.geometry.base not equal on grounded ice" 31 31 32 if any(md.geometry.base[n umpy.nonzero(md.mask.groundedice_levelset<=0.)]<md.geometry.bed[numpy.nonzero(md.mask.groundedice_levelset<=0.)]):32 if any(md.geometry.base[np.nonzero(md.mask.groundedice_levelset<=0.)]<md.geometry.bed[np.nonzero(md.mask.groundedice_levelset<=0.)]): 33 33 print "WARNING: md.geometry.base < md.geometry.bed on floating ice" 34 34 -
issm/trunk-jpl/src/m/parameterization/setmask.py
r20910 r21303 1 import numpy 1 import numpy as np 2 2 import os 3 3 from model import model … … 23 23 """ 24 24 #some checks on list of arguments 25 print type(md) 25 26 if not isinstance(md,model): 26 27 raise TypeError("setmask error message") … … 44 45 #arrays come from domain outlines that can intersect one another: 45 46 46 elementonfloatingice = n umpy.logical_and(elementonfloatingice,numpy.logical_not(elementongroundedice))47 elementongroundedice = n umpy.logical_not(elementonfloatingice)47 elementonfloatingice = np.logical_and(elementonfloatingice,np.logical_not(elementongroundedice)) 48 elementongroundedice = np.logical_not(elementonfloatingice) 48 49 49 50 #the order here is important. we choose vertexongroundedice as default on the grounding line. 50 vertexonfloatingice = n umpy.zeros(md.mesh.numberofvertices,'bool')51 vertexongroundedice = n umpy.zeros(md.mesh.numberofvertices,'bool')52 vertexongroundedice[md.mesh.elements[n umpy.nonzero(elementongroundedice),:]-1]=True53 vertexonfloatingice[n umpy.nonzero(numpy.logical_not(vertexongroundedice))]=True51 vertexonfloatingice = np.zeros(md.mesh.numberofvertices,'bool') 52 vertexongroundedice = np.zeros(md.mesh.numberofvertices,'bool') 53 vertexongroundedice[md.mesh.elements[np.nonzero(elementongroundedice),:]-1]=True 54 vertexonfloatingice[np.nonzero(np.logical_not(vertexongroundedice))]=True 54 55 #}}} 55 56 56 57 #level sets 57 md.mask.groundedice_levelset = -1.*n umpy.ones(md.mesh.numberofvertices)58 md.mask.groundedice_levelset[md.mesh.elements[n umpy.nonzero(elementongroundedice),:]-1]=1.58 md.mask.groundedice_levelset = -1.*np.ones(md.mesh.numberofvertices) 59 md.mask.groundedice_levelset[md.mesh.elements[np.nonzero(elementongroundedice),:]-1]=1. 59 60 60 61 if(len(args)): 61 md.mask.ice_levelset = 1.*n umpy.ones(md.mesh.numberofvertices)62 md.mask.ice_levelset = 1.*np.ones(md.mesh.numberofvertices) 62 63 icedomainfile = options.getfieldvalue('icedomain','none') 63 64 if not os.path.exists(icedomainfile): … … 65 66 #use contourtomesh to set ice values inside ice domain 66 67 vertexinsideicedomain,elementinsideicedomain=ContourToMesh(elements,x,y,icedomainfile,'node',1) 67 md.mask.ice_levelset[n umpy.nonzero(vertexinsideicedomain)[0]] = -1.68 md.mask.ice_levelset[np.nonzero(vertexinsideicedomain)[0]] = -1. 68 69 else: 69 md.mask.ice_levelset = -1.*n umpy.ones(md.mesh.numberofvertices)70 md.mask.ice_levelset = -1.*np.ones(md.mesh.numberofvertices) 70 71 71 72 return md -
issm/trunk-jpl/src/m/plot/applyoptions.py
r21297 r21303 1 import numpy as np1 import numpy as np 2 2 from cmaptools import truncate_colormap 3 3 from plot_contour import plot_contour -
issm/trunk-jpl/src/m/plot/checkplotoptions.py
r21283 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def checkplotoptions(md,options): -
issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 try: -
issm/trunk-jpl/src/m/plot/export_gl.py
r19056 r21303 2 2 from checkplotoptions import checkplotoptions 3 3 from model import model 4 import numpy as np4 import numpy as np 5 5 import math 6 6 from writejsfile import writejsfile -
issm/trunk-jpl/src/m/plot/plot_BC.py
r21283 r21303 4 4 print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled" 5 5 6 import numpy as np6 import numpy as np 7 7 from processmesh import processmesh 8 8 from applyoptions import applyoptions -
issm/trunk-jpl/src/m/plot/plot_icefront.py
r21283 r21303 3 3 except ImportError: 4 4 print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled" 5 import numpy as np5 import numpy as np 6 6 from processmesh import processmesh 7 7 from applyoptions import applyoptions -
issm/trunk-jpl/src/m/plot/plot_overlay.py
r21289 r21303 1 import numpy as np1 import numpy as np 2 2 from processmesh import processmesh 3 3 from processdata import processdata -
issm/trunk-jpl/src/m/plot/plot_quiver.py
r21283 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def plot_quiver(x,y,data,options,ax): -
issm/trunk-jpl/src/m/plot/plot_streamlines.py
r21254 r21303 1 import numpy as np1 import numpy as np 2 2 from processmesh import processmesh 3 3 from processdata import processdata -
issm/trunk-jpl/src/m/plot/plot_unit.py
r21288 r21303 5 5 import matplotlib as mpl 6 6 import matplotlib.pyplot as plt 7 import numpy as np7 import numpy as np 8 8 except ImportError: 9 9 print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled" -
issm/trunk-jpl/src/m/plot/plotmodel.py
r21298 r21303 1 import numpy as np1 import numpy as np 2 2 from plotoptions import plotoptions 3 3 from plotdoc import plotdoc -
issm/trunk-jpl/src/m/plot/processdata.py
r21283 r21303 1 import numpy as np1 import numpy as np 2 2 3 3 def processdata(md,data,options): -
issm/trunk-jpl/src/m/plot/processmesh.py
r21254 r21303 1 1 from math import isnan 2 import numpy as np2 import numpy as np 3 3 4 4 def processmesh(md,data,options): -
issm/trunk-jpl/src/m/plot/writejsfield.py
r19056 r21303 1 import numpy 1 import numpy as np 2 2 def writejsfield(fid,name,variable,nods): 3 3 #WRITEJSFIELD - write variable to javascript file … … 8 8 #write array: 9 9 #if not isinstance(variable, list): 10 if type(variable[0])==n umpy.float64:10 if type(variable[0])==np.float64: 11 11 fid.write('<!-- {0}{{{{{{-->\n'.format(name)) 12 12 fid.write('{0}=['.format(name)) -
issm/trunk-jpl/src/m/plot/writejsfile.py
r19056 r21303 1 import numpy 1 import numpy as np 2 2 from writejsfield import writejsfield 3 3 def writejsfile(filename,model,keyname): … … 48 48 fid.write('result["shortlabel"]="{0}";\n'.format(results[i].shortlabel)) 49 49 fid.write('result["unit"]="{0}";\n'.format(results[i].unit)) 50 if type(results[i].data)==n umpy.float64:50 if type(results[i].data)==np.float64: 51 51 fid.write('result["time_range"]=[{0},{1}];\n'.format(results[i].time_range[0],results[i].time_range[1])) 52 52 fid.write('results["{0}"]=result;\n'.format(i)) -
issm/trunk-jpl/src/m/solve/WriteData.py
r21144 r21303 1 import numpy 2 import math 1 import numpy as np 3 2 import struct 4 3 import pairoptions … … 44 43 if options.exist('scale'): 45 44 scale = options.getfieldvalue('scale') 46 if n umpy.size(data) > 1 :47 if n umpy.size(data,0)==timeserieslength:48 data=n umpy.array(data)45 if np.size(data) > 1 : 46 if np.size(data,0)==timeserieslength: 47 data=np.array(data) 49 48 data[0:-1,:] = scale*data[0:-1,:] 50 49 else: … … 52 51 else: 53 52 data = scale*data 54 if n umpy.size(data) > 1 :55 if n umpy.size(data,0)==timeserieslength:53 if np.size(data) > 1 : 54 if np.size(data,0)==timeserieslength: 56 55 yts = options.getfieldvalue('yts') 57 56 data[-1,:] = yts*data[-1,:] … … 119 118 120 119 if isinstance(data,bool): 121 data=n umpy.array([data])120 data=np.array([data]) 122 121 elif isinstance(data,(list,tuple)): 123 data=n umpy.array(data).reshape(-1,1)124 if n umpy.ndim(data) == 1:125 if n umpy.size(data):126 data=data.reshape(n umpy.size(data),1)122 data=np.array(data).reshape(-1,) 123 if np.ndim(data) == 1: 124 if np.size(data): 125 data=data.reshape(np.size(data),) 127 126 else: 128 127 data=data.reshape(0,0) … … 131 130 s=data.shape 132 131 #if matrix = NaN, then do not write anything 133 if s[0]==1 and s[1]==1 and math.isnan(data[0][0]):132 if np.ndim(data)==2 and np.product(s)==1 and np.all(np.isnan(data)): 134 133 s=(0,0) 135 134 136 135 #first write length of record 137 fid.write(struct.pack('i',4+4+8* s[0]*s[1]+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type136 fid.write(struct.pack('i',4+4+8*np.product(s)+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type 138 137 139 138 #write data code and matrix type: … … 142 141 143 142 #now write matrix 144 fid.write(struct.pack('i',s[0])) 145 fid.write(struct.pack('i',s[1])) 146 for i in xrange(s[0]): 147 for j in xrange(s[1]): 148 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 143 if np.ndim(data)==1: 144 fid.write(struct.pack('i',s[0])) 145 fid.write(struct.pack('i',1)) 146 for i in xrange(s[0]): 147 fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose 148 else: 149 fid.write(struct.pack('i',s[0])) 150 fid.write(struct.pack('i',s[1])) 151 for i in xrange(s[0]): 152 for j in xrange(s[1]): 153 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 149 154 # }}} 150 155 … … 152 157 153 158 if isinstance(data,(int,long)): 154 data=n umpy.array([data])159 data=np.array([data]) 155 160 elif isinstance(data,(list,tuple)): 156 data=n umpy.array(data).reshape(-1,1)157 if n umpy.ndim(data) == 1:158 if n umpy.size(data):159 data=data.reshape(n umpy.size(data),1)161 data=np.array(data).reshape(-1,) 162 if np.ndim(data) == 1: 163 if np.size(data): 164 data=data.reshape(np.size(data),) 160 165 else: 161 166 data=data.reshape(0,0) … … 164 169 s=data.shape 165 170 #if matrix = NaN, then do not write anything 166 if s[0]==1 and s[1]==1 and math.isnan(data[0][0]):171 if np.ndim(data)==2 and np.product(s)==1 and np.all(np.isnan(data)): 167 172 s=(0,0) 168 173 169 174 #first write length of record 170 fid.write(struct.pack('i',4+4+8* s[0]*s[1]+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type175 fid.write(struct.pack('i',4+4+8*np.product(s)+4+4)) #2 integers (32 bits) + the double matrix + code + matrix type 171 176 172 177 #write data code and matrix type: … … 175 180 176 181 #now write matrix 177 fid.write(struct.pack('i',s[0])) 178 fid.write(struct.pack('i',s[1])) 179 for i in xrange(s[0]): 180 for j in xrange(s[1]): 181 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 182 if np.ndim(data) == 1: 183 fid.write(struct.pack('i',s[0])) 184 fid.write(struct.pack('i',1)) 185 for i in xrange(s[0]): 186 fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose 187 else: 188 fid.write(struct.pack('i',s[0])) 189 fid.write(struct.pack('i',s[1])) 190 for i in xrange(s[0]): 191 for j in xrange(s[1]): 192 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 182 193 # }}} 183 194 … … 185 196 186 197 if isinstance(data,(bool,int,long,float)): 187 data=n umpy.array([data])198 data=np.array([data]) 188 199 elif isinstance(data,(list,tuple)): 189 data=n umpy.array(data).reshape(-1,1)190 if n umpy.ndim(data) == 1:191 if n umpy.size(data):192 data=data.reshape(n umpy.size(data),1)200 data=np.array(data).reshape(-1,) 201 if np.ndim(data) == 1: 202 if np.size(data): 203 data=data.reshape(np.size(data),) 193 204 else: 194 205 data=data.reshape(0,0) … … 197 208 s=data.shape 198 209 #if matrix = NaN, then do not write anything 199 if s[0]==1 and s[1]==1 and math.isnan(data[0][0]):210 if np.ndim(data)==1 and np.product(s)==1 and np.all(np.isnan(data)): 200 211 s=(0,0) 201 212 202 213 #first write length of record 203 recordlength=4+4+8* s[0]*s[1]+4+4; #2 integers (32 bits) + the double matrix + code + matrix type214 recordlength=4+4+8*np.product(s)+4+4; #2 integers (32 bits) + the double matrix + code + matrix type 204 215 if recordlength > 2**31 : 205 216 raise ValueError('field %s cannot be marshalled because it is larger than 4^31 bytes!' % enum) … … 212 223 213 224 #now write matrix 214 fid.write(struct.pack('i',s[0])) 215 fid.write(struct.pack('i',s[1])) 216 for i in xrange(s[0]): 217 for j in xrange(s[1]): 218 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 225 if np.ndim(data) == 1: 226 fid.write(struct.pack('i',s[0])) 227 fid.write(struct.pack('i',1)) 228 for i in xrange(s[0]): 229 fid.write(struct.pack('d',float(data[i]))) #get to the "c" convention, hence the transpose 230 else: 231 fid.write(struct.pack('i',s[0])) 232 fid.write(struct.pack('i',s[1])) 233 for i in xrange(s[0]): 234 for j in xrange(s[1]): 235 fid.write(struct.pack('d',float(data[i][j]))) #get to the "c" convention, hence the transpose 219 236 # }}} 220 237 … … 225 242 for matrix in data: 226 243 if isinstance(matrix,(bool,int,long,float)): 227 matrix=n umpy.array([matrix])244 matrix=np.array([matrix]) 228 245 elif isinstance(matrix,(list,tuple)): 229 matrix=n umpy.array(matrix).reshape(-1,1)230 if n umpy.ndim(matrix) == 1:231 if n umpy.size(matrix):232 matrix=matrix.reshape(n umpy.size(matrix),1)246 matrix=np.array(matrix).reshape(-1,) 247 if np.ndim(matrix) == 1: 248 if np.size(matrix): 249 matrix=matrix.reshape(np.size(matrix),) 233 250 else: 234 251 matrix=matrix.reshape(0,0) 235 252 236 253 s=matrix.shape 237 recordlength+=4*2+ s[0]*s[1]*8 #row and col of matrix + matrix of doubles254 recordlength+=4*2+np.product(s)*8 #row and col of matrix + matrix of doubles 238 255 239 256 #write length of record … … 249 266 for matrix in data: 250 267 if isinstance(matrix,(bool,int,long,float)): 251 matrix=n umpy.array([matrix])268 matrix=np.array([matrix]) 252 269 elif isinstance(matrix,(list,tuple)): 253 matrix=n umpy.array(matrix).reshape(-1,1)254 if n umpy.ndim(matrix) == 1:255 matrix=matrix.reshape(n umpy.size(matrix),1)270 matrix=np.array(matrix).reshape(-1,) 271 if np.ndim(matrix) == 1: 272 matrix=matrix.reshape(np.size(matrix),) 256 273 257 274 s=matrix.shape 258 fid.write(struct.pack('i',s[0])) 259 fid.write(struct.pack('i',s[1])) 260 for i in xrange(s[0]): 261 for j in xrange(s[1]): 262 fid.write(struct.pack('d',float(matrix[i][j]))) 275 if np.ndim(data) == 1: 276 fid.write(struct.pack('i',s[0])) 277 fid.write(struct.pack('i',1)) 278 for i in xrange(s[0]): 279 fid.write(struct.pack('d',float(matrix[i]))) #get to the "c" convention, hence the transpose 280 else: 281 fid.write(struct.pack('i',s[0])) 282 fid.write(struct.pack('i',s[1])) 283 for i in xrange(s[0]): 284 for j in xrange(s[1]): 285 fid.write(struct.pack('d',float(matrix[i][j]))) 263 286 # }}} 264 287 -
issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
r21267 r21303 1 1 import struct 2 import numpy 2 import numpy as np 3 3 from collections import OrderedDict 4 4 import results as resultsclass … … 153 153 M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] 154 154 if type==1: 155 field=n umpy.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float)155 field=np.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float) 156 156 elif type==2: 157 157 field=struct.unpack('%ds' % M,fid.read(M))[0][:-1] … … 159 159 N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] 160 160 # field=transpose(fread(fid,[N M],'double')); 161 field=n umpy.zeros(shape=(M,N),dtype=float)161 field=np.zeros(shape=(M,N),dtype=float) 162 162 for i in xrange(M): 163 163 field[i,:]=struct.unpack('%dd' % N,fid.read(N*struct.calcsize('d'))) -
issm/trunk-jpl/src/m/solve/solve.py
r21097 r21303 7 7 from waitonlock import waitonlock 8 8 from loadresultsfromcluster import loadresultsfromcluster 9 import MatlabFuncs as m10 9 11 10 def solve(md,solutionstring,*args): … … 76 75 md.private.solution=solutionstring 77 76 cluster=md.cluster 78 if m.strcmpi(options.getfieldvalue('batch','no'),'yes'):77 if options.getfieldvalue('batch','no')=='yes': 79 78 batch=1 80 79 else: … … 82 81 83 82 #check model consistency 84 if m.strcmpi(options.getfieldvalue('checkconsistency','yes'),'yes'):83 if options.getfieldvalue('checkconsistency','yes')=='yes': 85 84 print "checking model consistency" 86 85 ismodelselfconsistent(md) … … 117 116 118 117 #Stop here if batch mode 119 if m.strcmpi(options.getfieldvalue('batch','no'),'yes'):118 if options.getfieldvalue('batch','no')=='yes': 120 119 print 'batch mode requested: not launching job interactively' 121 120 print 'launch solution sequence on remote cluster by hand'
Note:
See TracChangeset
for help on using the changeset viewer.