Changeset 13470

Timestamp:

09/27/12 14:46:49 (12 years ago)

Author:

jschierm

Message:

CHG: More changes relating to python 1-d vectors vs. 2-d arrays (plus some cleanup).

Location:

issm/trunk-jpl/src/m

Files:

• boundaryconditions/SetIceShelfBC.py (modified) (3 diffs)
• boundaryconditions/SetMarineIceSheetBC.m (modified) (2 diffs)
• boundaryconditions/SetMarineIceSheetBC.py (added)
• classes/model/model.py (modified) (1 diff)
• geometry/FlagElements.py (modified) (2 diffs)
• mesh/bamg.py (modified) (2 diffs)

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

@@ -26 +26 @@
                         raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
                 [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
-                nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront).astype(float)
+                nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
         else:
-                nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1))
+                nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
 
 #       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
-        pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif]
-        md.diagnostic.spcvx=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
-        md.diagnostic.spcvy=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
-        md.diagnostic.spcvz=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]
+        md.diagnostic.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.diagnostic.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.diagnostic.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
         md.diagnostic.spcvx[pos]=0
         md.diagnostic.spcvy[pos]=0
         md.diagnostic.spcvz[pos]=0
-        md.diagnostic.referential=float('NaN')*numpy.ones((md.mesh.numberofvertices,6))
+        md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
 
         #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:
-                print '      boundary conditions for diagnostic model: spc set as observed velocities'
+                print "      boundary conditions for diagnostic model: spc set as observed velocities"
                 md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
                 md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos]
         else:
-                print '      boundary conditions for diagnostic model: spc set as zero'
+                print "      boundary conditions for diagnostic model: spc set as zero"
 
         #segment on Ice Front
         #segment on Neumann (Ice Front)
 #       pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
-        pos=[i for i,(noif1,noif2) in enumerate(zip(nodeonicefront[md.mesh.segments[:,0].astype('int')-1],nodeonicefront[md.mesh.segments[:,1].astype('int')-1])) if noif1 or noif2]
+        pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
         if   md.mesh.dimension==2:
                 pressureload=md.mesh.segments[pos,:]
@@ -60 +60 @@
                 for i in xrange(1,md.mesh.numberoflayers):
 #                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
-                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
+                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
 
         #Add water or air enum depending on the element
 #       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
-        pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))))
+        pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
 
         #plug onto model
@@ -72 +72 @@
         if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
                 md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-                print '      no surfaceforcings.precipitation specified: values set as zero'
+                print "      no surfaceforcings.precipitation specified: values set as zero"
         if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
                 md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
-                print '      no surfaceforcings.mass_balance specified: values set as zero'
+                print "      no surfaceforcings.mass_balance specified: values set as zero"
         if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
                 md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-                print '      no basalforcings.melting_rate specified: values set as zero'
+                print "      no basalforcings.melting_rate specified: values set as zero"
         if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
                 md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
-                print '      no balancethickness.thickening_rate specified: values set as zero'
+                print "      no balancethickness.thickening_rate specified: values set as zero"
 
-        md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
-        md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 
-        if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices:
-                md.thermal.spctemperature=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        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))
 #               pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-                pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos]
-                md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    # impose observed temperature on surface
-                if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
+                pos=numpy.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))
         else:
-                print '      no thermal boundary conditions created: no observed temperature found'
+                print "      no thermal boundary conditions created: no observed temperature found"
 
         return md

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

@@ -53 +53 @@
 
 md.hydrology.spcwatercolumn=zeros(md.mesh.numberofvertices,2);
-pos=find(md.mesh.vertexonboundary); 
+pos=find(md.mesh.vertexonboundary);
 md.hydrology.spcwatercolumn(pos,1)=1;
 
@@ -101 +101 @@
         if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
                 md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
-                md.basalforcings.geothermalflux(find(md.mask.vertexongroundedice))=50*10^-3; %50mW/m2
+                md.basalforcings.geothermalflux(find(md.mask.vertexongroundedice))=50.*10.^-3; %50mW/m2
         end
 else

issm/trunk-jpl/src/m/classes/model/model.py

@@ -167 +167 @@
                 print ("model not consistent: %s" % string)
                 self.private.isconsistent=False
+                return self
         # }}}
 

issm/trunk-jpl/src/m/geometry/FlagElements.py

@@ -41 +41 @@
                                 if len(region)>3 and not strcmp(region[-4:],'.exp'):
                                         raise IOError("Error: File 'region' not found!" % region)
-                                raise RuntimeError("FlagElements -- basinzoom not yet converted.")
+                                raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.")
                                 xlim,ylim=basinzoom('basin',region)
                                 flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0])).astype(float)
@@ -47 +47 @@
                         else:
                                 #ok, flag elements
-                                [flag,fnone]=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1)
+                                [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1)
 
                 if invert:

issm/trunk-jpl/src/m/mesh/bamg.py

@@ -121 +121 @@
 
                                 elif numpy.any(numpy.logical_not(flags)):
-                                        raise RuntimeError("bamg.m for rifts is not complete.")
+                                        raise RuntimeError("bamg.py for rifts is not complete.")
                                         #We LOTS of work to do
                                         print("Rift tip outside of or on the domain has been detected and is being processed...")
@@ -347 +347 @@
 def processgeometry(geom,tol,outline):    # {{{
 
-        raise RuntimeError("bamg.m/processgeometry is not complete.")
+        raise RuntimeError("bamg.py/processgeometry is not complete.")
         #Deal with edges
         print("Checking Edge crossing...")