Changeset 13975 for issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py

Timestamp:

Nov 16, 2012, 8:10:16 AM (14 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 13974

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/boundaryconditions/SetIceShelfBC.py (modified) (4 diffs)

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

@@ -1 +1 @@
 import os
 import numpy
+from ContourToMesh import *
 
 def SetIceShelfBC(md,icefrontfile=''):
@@ -6 +7 @@
         SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
 
-           Neumann BC are used on the ice front (an ANRGUS contour around the ice front
+           Neumann BC are used on the ice front (an ARGUS contour around the ice front
            must be given in input)
            Dirichlet BC are used elsewhere for diagnostic
@@ -24 +25 @@
                 if not os.path.exists(icefrontfile):
                         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=double(md.mesh.vertexonboundary and nodeinsideicefront)
+                [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.reshape(-1)).astype(float)
         else:
-                nodeonicefront=numpy.zeros(md.mesh.numberofvertices)
+                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)
-        md.diagnostic.spcvy=float('NaN')*numpy.ones(md.mesh.numberofvertices)
-        md.diagnostic.spcvz=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+        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 numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices:
-                print '      boundary conditions for diagnostic model: spc set as observed velocities'
+        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:
+                #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)
+                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,:]
         elif md.mesh.dimension==3:
 #               pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
-                pressureload_layer1=numpy.concatenate((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]),axis=1)
+                pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
                 pressureload=numpy.zeros((0,5))
                 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.concatenate((pressureload,numpy.concatenate((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d),axis=1)),axis=0)
+                        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.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)
+        pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
 
         #plug onto model
@@ -70 +76 @@
         #Create zeros basalforcings and surfaceforcings
         if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
-                md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices)
-                print '      no surfaceforcings.precipitation specified: values set as zero'
+                md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
+                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)
-                print '      no surfaceforcings.mass_balance specified: values set as zero'
+                md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
+                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)
-                print '      no basalforcings.melting_rate specified: values set as zero'
+                md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
+                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)
-                print '      no balancethickness.thickening_rate specified: values set as zero'
+                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+                print "      no balancethickness.thickening_rate specified: values set as zero"
 
-        md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
-        md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+        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)
+        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:
-                        md.basalforcings.geothermalflux=numpy.zeros(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