source: issm/oecreview/Archive/15392-16133/ISSM-15994-15995.diff@ 16134

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

@@ -26 +26 @@
                 #User provided Front.exp, use it
                 if not os.path.exists(icefrontfile):
                         raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
-                [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
+                [incontour,dum]=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
+                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))
         else:
                 #Guess where the ice front is
-                vertexonfloatingice=numpy.logical(md.mask.groundedice_levelset<0.)
-                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice)
+                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);
         pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]