Changeset 21303 for issm/trunk-jpl/src/m/boundaryconditions

Timestamp:

10/21/16 08:49:08 (9 years ago)

Author:

bdef

Message:

CHG: major cosmetic clean up we now ask shape (N,) rather than (N,1) also uniformised numpy imports

Location:

issm/trunk-jpl/src/m/boundaryconditions

Files:

• PattynSMB.py (modified) (1 diff)
• SetIceSheetBC.py (modified) (3 diffs)
• SetIceShelfBC.py (modified) (5 diffs)
• SetMarineIceSheetBC.py (modified) (6 diffs)
• love_numbers.py (modified) (1 diff)

issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py

@@ -1 +1 @@
 import os
-import numpy as np
+import numpy as  np
 def PattynSMB(md,Tf):
         """

issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py

@@ -1 +1 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 
@@ -14 +14 @@
 
         #node on Dirichlet
-        pos=numpy.nonzero(md.mesh.vertexonboundary)
-        md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
+        pos=np.nonzero(md.mesh.vertexonboundary)
+        md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
         md.stressbalance.spcvx[pos]=0
         md.stressbalance.spcvy[pos]=0
         md.stressbalance.spcvz[pos]=0
-        md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-        md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+        md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+        md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
         #Dirichlet Values
-        if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+        if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
                 print "      boundary conditions for stressbalance model: spc set as observed velocities"
                 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -39 +39 @@
 
         #Deal with other boundary conditions
-        if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+        if np.all(np.isnan(md.balancethickness.thickening_rate)):
+                md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
                 print "      no balancethickness.thickening_rate specified: values set as zero"
-        md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-        if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-                md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+                md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
                 if hasattr(md.mesh,'vertexonsurface'):
-                        pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+                        pos=np.nonzero(md.mesh.vertexonsurface)[0]
                         md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-                if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
-                        md.basalforcings.geothermalflux=50.*10**-3*numpy.ones((md.mesh.numberofvertices,1))    #50 mW/m^2
+                if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
+                        md.basalforcings.geothermalflux=50.*10**-3*np.ones((md.mesh.numberofvertices))    #50 mW/m^2
         else:
                 print "      no thermal boundary conditions created: no observed temperature found"

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

@@ -1 +1 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 import MatlabFuncs as m
@@ -27 +27 @@
                         raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
                 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-                nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
+                nodeonicefront=np.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
         else:
-                nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
+                nodeonicefront=np.zeros((md.mesh.numberofvertices),bool)
 
 #       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
-        pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]
-        md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-        md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+        pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(nodeonicefront)))[0]
+        md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+        md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
         #Icefront position
-        pos=numpy.nonzero(nodeonicefront)[0]
+        pos=np.nonzero(nodeonicefront)[0]
         md.mask.ice_levelset[pos]=0
 
@@ -53 +53 @@
                 values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
                 segmentsfront=1-values
-                numpy.sum(segmentsfront,axis=1)!=numbernodesfront
-                segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
+                np.sum(segmentsfront,axis=1)!=numbernodesfront
+                segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
                 #Find all nodes for these segments and spc them
                 pos=md.mesh.segments[segments,0:-1]-1
         else:
-                pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+                pos=np.nonzero(md.mesh.vertexonboundary)[0]
         md.stressbalance.spcvx[pos]=0
         md.stressbalance.spcvy[pos]=0
@@ -64 +64 @@
                                                                                                                                                                                                                                            
         #Dirichlet Values
-        if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+        if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
                 #reshape to rank-2 if necessary to match spc arrays
-                if numpy.ndim(md.inversion.vx_obs)==1:
-                        md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,1)
-                if numpy.ndim(md.inversion.vy_obs)==1:
-                        md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,1)
+                if np.ndim(md.inversion.vx_obs)==1:
+                        md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,)
+                if np.ndim(md.inversion.vy_obs)==1:
+                        md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,)
                 print "      boundary conditions for stressbalance model: spc set as observed velocities"
                 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -81 +81 @@
 
         #Deal with other boundary conditions
-        if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+        if np.all(np.isnan(md.balancethickness.thickening_rate)):
+                md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
                 print "      no balancethickness.thickening_rate specified: values set as zero"
-        md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-        if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-                md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+                md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
                 if hasattr(md.mesh,'vertexonsurface'):
-                        pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+                        pos=np.nonzero(md.mesh.vertexonsurface)[0]
                         md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-                if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
-                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
+                if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
+                        md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices))
         else:
                 print "      no thermal boundary conditions created: no observed temperature found"

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

@@ -1 +1 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 import MatlabFuncs as m
@@ -29 +29 @@
                         raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
                 incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))
+                vertexonicefront=np.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))
         else:
                 #Guess where the ice front is
-                vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,1))
-                pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]
+                vertexonfloatingice=np.zeros((md.mesh.numberofvertices))
+                pos=np.nonzero(np.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]
                 vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1.
-                vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.)
+                vertexonicefront=np.logical_and(np.reshape(md.mesh.vertexonboundary,(-1,)),vertexonfloatingice>0.)
 
 #       pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
-        pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]
-        if not numpy.size(pos):
+        pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(vertexonicefront)))[0]
+        if not np.size(pos):
                 print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually."
 
-        md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
-        md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-        md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+        md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
+        md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+        md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
         #Position of ice front
-        pos=numpy.nonzero(vertexonicefront)[0]
+        pos=np.nonzero(vertexonicefront)[0]
         md.mask.ice_levelset[pos]=0
 
@@ -62 +62 @@
                 values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
                 segmentsfront=1-values
-                numpy.sum(segmentsfront,axis=1)!=numbernodesfront
-                segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
+                np.sum(segmentsfront,axis=1)!=numbernodesfront
+                segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
                 #Find all nodes for these segments and spc them
                 pos=md.mesh.segments[segments,0:-1]-1
         else:
-                pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+                pos=np.nonzero(md.mesh.vertexonboundary)[0]
         md.stressbalance.spcvx[pos]=0
         md.stressbalance.spcvy[pos]=0
@@ -73 +73 @@
 
         #Dirichlet Values
-        if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+        if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
                 print "      boundary conditions for stressbalance model: spc set as observed velocities"
                 md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -80 +80 @@
                 print "      boundary conditions for stressbalance model: spc set as zero"
 
-        md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2))
-        pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+        md.hydrology.spcwatercolumn=np.zeros((md.mesh.numberofvertices,2))
+        pos=np.nonzero(md.mesh.vertexonboundary)[0]
         md.hydrology.spcwatercolumn[pos,0]=1
 
@@ -89 +89 @@
 
         #Deal with other boundary conditions
-        if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+        if np.all(np.isnan(md.balancethickness.thickening_rate)):
+                md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
                 print "      no balancethickness.thickening_rate specified: values set as zero"
 
-        md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-        md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+        md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-        if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-                md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+                md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
                 if hasattr(md.mesh,'vertexonsurface'):
-                        pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+                        pos=np.nonzero(md.mesh.vertexonsurface)[0]
                         md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-                if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
-                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
-                        md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
+                if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
+                        md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices))
+                        md.basalforcings.geothermalflux[np.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
         else:
                 print "      no thermal boundary conditions created: no observed temperature found"

issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py

@@ -1 +1 @@
 from MatlabFuncs import *
 from model import *
-from numpy import *
+from np.import *
 
 def love_numbers(value,*varargin):