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

Index: ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py        (revision 15994)
+++ ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py        (revision 15995)
@@ -26,12 +26,14 @@
                #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]