Changeset 12191

Timestamp:

05/03/12 12:43:52 (13 years ago)

Author:

Eric.Larour

Message:

Further translation of SetIceShelfBC.m to python

File:

• issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py (modified) (3 diffs)

issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py

@@ -1 @@
-def SetIceShelfBC(md,icefrontfile='')
+from numpy import *
+def SetIceShelfBC(md,icefrontfile=''):
         #SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
         #
@@ -22 +23 @@
                         print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found'
                         return []
-                nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2);
-                nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront);
+                nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
+                nodeonicefront=double(md.mesh.vertexonboundary.astype(bool) and nodeinsideicefront.astype(bool))
         
-        pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
-        md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
-        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
-        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-        md.diagnostic.spcvx(pos)=0;
-        md.diagnostic.spcvy(pos)=0;
-        md.diagnostic.spcvz(pos)=0;
-        md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
+        pos=argwhere(logical_and(md.mesh.vertexonboundary.astype(bool),~nodeonicefront.astype(bool)))
+        md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices);
+        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices);
+        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices);
+        md.diagnostic.spcvx[pos]=0;
+        md.diagnostic.spcvy[pos]=0;
+        md.diagnostic.spcvz[pos]=0;
+        md.diagnostic.referential=float('NaN')*ones((md.mesh.numberofvertices,6),float);
 
         #Dirichlet Values
-        if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
-                disp('      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
-                disp('      boundary conditions for diagnostic model: spc set as zero');
-        end
+        if ~isnan(md.inversion.vx_obs) and ~isnan(md.inversion.vy_obs):
+                if (len(md.inversion.vx_obs)==md.mesh.numberofvertices) and (len(md.inversion.vy_obs)==md.mesh.numberofvertices):
+                        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'
 
         #segment on Ice Front
         #segment on Neumann (Ice Front)
-        pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
-        if (md.mesh.dimension==2)
-                pressureload=md.mesh.segments(pos,:);
-        elseif 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)];
+        segs1=md.mesh.segments[:,0].astype(int)-1
+        segs2=md.mesh.segments[:,1].astype(int)-1
+
+        pos=argwhere(logical_and(nodeonicefront[segs1].astype(bool),nodeonicefront[segs2].astype(bool)))
+        if (md.mesh.dimension==2):
+                pressureload=md.mesh.segments[pos,:];
+        elif md.mesh.dimension==3:
+                pressureload_layer=[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=[];
-                for i=1:md.mesh.numberoflayers-1,
-                        pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
-                end
-        end
+                for i in range(1,md.mesh.numberoflayers-1):
+                        pressureload=[pressureload ,pressureload_layer1[:,1:4]+(i-1)*md.mesh.numberofvertices2d, pressureload_layer1[:,5]+(i-1)*md.mesh.numberofelements2d ];
 
         #Add water or air enum depending on the element
-        pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
+        pressureload=[pressureload, 1*md.mask.elementonfloatingice(pressureload[:,end])];
 
         #plug onto model
@@ -63 +65 @@
 
         #Create zeros basalforcings and surfaceforcings
-        if isnan(md.surfaceforcings.precipitation),
-                md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
-                disp('      no surfaceforcings.precipitation specified: values set as zero');
-        end
-        if isnan(md.surfaceforcings.mass_balance),
-                md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-                disp('      no surfaceforcings.mass_balance specified: values set as zero');
-        end
-        if isnan(md.basalforcings.melting_rate),
-                md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-                disp('      no basalforcings.melting_rate specified: values set as zero');
-        end
-        if isnan(md.balancethickness.thickening_rate),
-                md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
-                disp('      no balancethickness.thickening_rate specified: values set as zero');
-        end
+        if isnan(md.surfaceforcings.precipitation):
+                md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices);
+                print '      no surfaceforcings.precipitation specified: values set as zero'
+        if isnan(md.surfaceforcings.mass_balance):
+                md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices);
+                print '      no surfaceforcings.mass_balance specified: values set as zero'
+        if isnan(md.basalforcings.melting_rate):
+                md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices)
+                print '      no basalforcings.melting_rate specified: values set as zero'
+        if isnan(md.balancethickness.thickening_rate):
+                md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices);
+                print '      no balancethickness.thickening_rate specified: values set as zero'
 
-        md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
-        md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
+        md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices);
+        md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices);
 
-        if (length(md.initialization.temperature)==md.mesh.numberofvertices),
-                md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
-                pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); #impose observed temperature on surface
-                if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
-                        md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
-                end
-        else
-                disp('      no thermal boundary conditions created: no observed temperature found');
-        end
+        if (len(md.initialization.temperature)==md.mesh.numberofvertices):
+                md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices);
+                pos=argwhere(md.mesh.vertexonsurface); md.thermal.spctemperature[pos]=md.initialization.temperature[pos]; #impose observed temperature on surface
+                if (len(md.basalforcings.geothermalflux) !=md.mesh.numberofvertices):
+                        md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices);
+        else:
+                print '      no thermal boundary conditions created: no observed temperature found'
+
+        return md