Changeset 17766

Timestamp:

04/18/14 09:32:00 (11 years ago)

Author:

cborstad

Message:

CHG: reverting change to postpone python handling of vectors/arrays for now

File:

• issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (modified) (4 diffs)

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

@@ -32 +32 @@
         else:
                 #Guess where the ice front is
-                vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,))
-                pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0,axis=1) > 0)[0]
-                #vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1.
-                vertexonfloatingice[md.mesh.elements[pos]-1]=1.
-                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice>0.)
+                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[md.mesh.elements[pos].astype(int)-1]=1.
+                vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.)
 
 #       pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
@@ -43 +42 @@
                 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.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+        md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
         md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
         md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
@@ -69 +68 @@
         else:
                 pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
-        md.stressbalance.spcvx[pos]=0
-        md.stressbalance.spcvy[pos]=0
-        md.stressbalance.spcvz[pos]=0
+        md.stressbalance.spcvx[pos[:]]=0
+        md.stressbalance.spcvy[pos[:]]=0
+        md.stressbalance.spcvz[pos[:]]=0
 
         #Dirichlet Values
@@ -91 +90 @@
         #Deal with other boundary conditions
         if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,))
+                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
                 print "      no balancethickness.thickening_rate specified: values set as zero"
 
-        md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-        md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,))
+        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))
 
         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,))
+                md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
                 if hasattr(md.mesh,'vertexonsurface'):
                         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,))
-                        md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=50.*10.**-3    #50mW/m2
+                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
+                        md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
         else:
                 print "      no thermal boundary conditions created: no observed temperature found"
 
         return md
+