Index: /issm/trunk-jpl/src/m/exp/contourlevelzero.py
===================================================================
--- /issm/trunk-jpl/src/m/exp/contourlevelzero.py	(revision 26374)
+++ /issm/trunk-jpl/src/m/exp/contourlevelzero.py	(revision 26374)
@@ -0,0 +1,237 @@
+import os.path
+import numpy as np
+from collections import OrderedDict
+
+def contourlevelzero(md,mask,level):
+    """CONTOURLEVELZERO - figure out the zero level (or offset thereof, specified by the level value)
+                       of a vectorial mask, and vectorialize it into an exp or shp compatible structure.
+    
+       Usage:
+          contours=contourlevelzero(md,mask,level)
+    
+       See also: PLOT_CONTOUR
+    """
+    
+    #process data 
+    if md.mesh.dimension()==3:
+        x = md.mesh.x2d
+        y = md.mesh.y2d
+        z=md.mesh.z
+        index=md.mesh.elements2d-1
+    else:
+        x=md.mesh.x
+        y=md.mesh.y
+        index=md.mesh.elements-1
+        z=np.zeros((md.mesh.numberofvertices,1))
+        
+    if len(mask)==0:
+        raise OSError("mask provided is empty")
+    
+    if md.mesh.dimension()==3:
+        if len(mask)!=md.mesh.numberofvertices2d: 
+            raise OSError("mask provided should be specified at the vertices of the mesh")
+    else:
+        if len(mask)!=md.mesh.numberofvertices:
+            raise OSError("mask provided should be specified at the vertices of the mesh")
+        
+    #initialization of some variables
+    numberofelements=np.size(index,0)
+    elementslist=np.c_[0:numberofelements]
+    c=[]
+    h=[]
+    
+    #get unique edges in mesh
+    #1: list of edges
+    edges=np.vstack((np.vstack((index[:,(0,1)],index[:,(1,2)])),index[:,(2,0)]))
+
+    #2: find unique edges
+    [edges,J]=np.unique(np.sort(edges,1),axis=0,return_inverse=True)
+    #3: unique edge numbers
+    vec=J
+    #4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
+    #   the same edge number)
+    edges_tria=np.hstack((np.hstack((vec[elementslist],vec[elementslist+numberofelements])),vec[elementslist+2*numberofelements]))
+    
+    #segments [nodes1 nodes2]
+    Seg1=index[:,(0,1)]
+    Seg2=index[:,(1,2)]
+    Seg3=index[:,(2,0)]
+    
+    #segment numbers [1;4;6;...]
+    Seg1_num=edges_tria[:,0]
+    Seg2_num=edges_tria[:,1]
+    Seg3_num=edges_tria[:,2]
+    
+    #value of data on each tips of the segments
+    Data1=mask[Seg1]
+    Data2=mask[Seg2]
+    Data3=mask[Seg3]
+    
+    #get the ranges for each segment
+    Range1=np.sort(Data1,1)
+    Range2=np.sort(Data2,1)
+    Range3=np.sort(Data3,1)
+    
+    #find the segments that contain this value
+    pos1=(Range1[:,0]<level) & (Range1[:,1]>=level)
+    pos2=(Range2[:,0]<level) & (Range2[:,1]>=level)
+    pos3=(Range3[:,0]<level) & (Range3[:,1]>=level)
+    
+    #get elements
+    poselem12=(pos1) & (pos2)
+    poselem13=(pos1) & (pos3)
+    poselem23=(pos2) & (pos3)
+    poselem=np.where((poselem12) | (poselem13) | (poselem23))
+    poselem=poselem[0]
+    numelems=len(poselem)
+    
+    #if no element has been flagged, skip to the next level
+    if numelems==0:
+        raise Exception('contourlevelzero warning message: no elements found with corresponding level value in mask')
+        contours=[]
+        return contours
+    
+    #go through the elements and build the coordinates for each segment (1 by element)
+    x1=np.zeros((numelems,1))
+    x2=np.zeros((numelems,1))
+    y1=np.zeros((numelems,1))
+    y2=np.zeros((numelems,1))
+    z1=np.zeros((numelems,1))
+    z2=np.zeros((numelems,1))
+    
+    edge_l=np.zeros((numelems,2))
+    
+    for j in range(0,numelems):
+        
+        with np.errstate(divide='ignore', invalid='ignore'):
+            weight1=np.divide(level-Data1[poselem[j],0],Data1[poselem[j],1]-Data1[poselem[j],0])
+            weight2=np.divide(level-Data2[poselem[j],0],Data2[poselem[j],1]-Data2[poselem[j],0])
+            weight3=np.divide(level-Data3[poselem[j],0],Data3[poselem[j],1]-Data3[poselem[j],0])
+        
+        if poselem12[poselem[j]]==True:
+            
+            x1[j]=x[Seg1[poselem[j],0]]+weight1*[x[Seg1[poselem[j],1]]-x[Seg1[poselem[j],0]]]
+            x2[j]=x[Seg2[poselem[j],0]]+weight2*[x[Seg2[poselem[j],1]]-x[Seg2[poselem[j],0]]]
+            y1[j]=y[Seg1[poselem[j],0]]+weight1*[y[Seg1[poselem[j],1]]-y[Seg1[poselem[j],0]]]
+            y2[j]=y[Seg2[poselem[j],0]]+weight2*[y[Seg2[poselem[j],1]]-y[Seg2[poselem[j],0]]]
+            z1[j]=z[Seg1[poselem[j],0]]+weight1*[z[Seg1[poselem[j],1]]-z[Seg1[poselem[j],0]]]
+            z2[j]=z[Seg2[poselem[j],0]]+weight2*[z[Seg2[poselem[j],1]]-z[Seg2[poselem[j],0]]]
+            
+            edge_l[j,0]=Seg1_num[poselem[j]]
+            edge_l[j,1]=Seg2_num[poselem[j]]
+        elif poselem13[poselem[j]]==True:
+            
+            x1[j]=x[Seg1[poselem[j],0]]+weight1*[x[Seg1[poselem[j],1]]-x[Seg1[poselem[j],0]]]
+            x2[j]=x[Seg3[poselem[j],0]]+weight3*[x[Seg3[poselem[j],1]]-x[Seg3[poselem[j],0]]]
+            y1[j]=y[Seg1[poselem[j],0]]+weight1*[y[Seg1[poselem[j],1]]-y[Seg1[poselem[j],0]]]
+            y2[j]=y[Seg3[poselem[j],0]]+weight3*[y[Seg3[poselem[j],1]]-y[Seg3[poselem[j],0]]]
+            z1[j]=z[Seg1[poselem[j],0]]+weight1*[z[Seg1[poselem[j],1]]-z[Seg1[poselem[j],0]]]
+            z2[j]=z[Seg3[poselem[j],0]]+weight3*[z[Seg3[poselem[j],1]]-z[Seg3[poselem[j],0]]]
+            
+            edge_l[j,0]=Seg1_num[poselem[j]]
+            edge_l[j,1]=Seg3_num[poselem[j]]
+        elif poselem23[poselem[j]]==True:
+            
+            x1[j]=x[Seg2[poselem[j],0]]+weight2*[x[Seg2[poselem[j],1]]-x[Seg2[poselem[j],0]]]
+            x2[j]=x[Seg3[poselem[j],0]]+weight3*[x[Seg3[poselem[j],1]]-x[Seg3[poselem[j],0]]]
+            y1[j]=y[Seg2[poselem[j],0]]+weight2*[y[Seg2[poselem[j],1]]-y[Seg2[poselem[j],0]]]
+            y2[j]=y[Seg3[poselem[j],0]]+weight3*[y[Seg3[poselem[j],1]]-y[Seg3[poselem[j],0]]]
+            z1[j]=z[Seg2[poselem[j],0]]+weight2*[z[Seg2[poselem[j],1]]-z[Seg2[poselem[j],0]]]
+            z2[j]=z[Seg3[poselem[j],0]]+weight3*[z[Seg3[poselem[j],1]]-z[Seg3[poselem[j],0]]]
+
+            edge_l[j,0]=Seg2_num[poselem[j]]
+            edge_l[j,1]=Seg3_num[poselem[j]]
+
+        #else:
+	    #it shoud not go here
+            
+    #now that we have the segments, we must try to connect them...
+    
+    #loop over the subcontours
+    contours=[]
+    
+    while len(edge_l)>0:
+        
+        #take the right edge of the second segment and connect it to the next segments if any
+        e1=edge_l[0,0]
+        e2=edge_l[0,1]
+        xc=np.vstack((x1[0],x2[0]))
+        yc=np.vstack((y1[0],y2[0]))
+        zc=np.vstack((z1[0],z2[0]))
+        #erase the lines corresponding to this edge
+        edge_l=np.delete(edge_l,0,axis=0)
+        x1=np.delete(x1,0,axis=0)
+        x2=np.delete(x2,0,axis=0)
+        y1=np.delete(y1,0,axis=0)
+        y2=np.delete(y2,0,axis=0)
+        z1=np.delete(z1,0,axis=0)
+        z2=np.delete(z2,0,axis=0)
+        pos1=np.where(edge_l==e1)
+        
+        while len(pos1[0])>0:
+            
+            if np.all(pos1[1]==0):
+                xc=np.vstack((x2[pos1[0]],xc))
+                yc=np.vstack((y2[pos1[0]],yc))
+                zc=np.vstack((z2[pos1[0]],zc))
+                #next edge:
+                e1=edge_l[pos1[0],1]
+            else:
+                xc=np.vstack((x1[pos1[0]],xc))
+                yc=np.vstack((y1[pos1[0]],yc))
+                zc=np.vstack((z1[pos1[0]],zc))
+                #next edge:
+                e1=edge_l[pos1[0],0]
+                
+            #erase the lines of this
+            edge_l=np.delete(edge_l,pos1[0],axis=0)
+            x1=np.delete(x1,pos1[0],axis=0)
+            x2=np.delete(x2,pos1[0],axis=0)
+            y1=np.delete(y1,pos1[0],axis=0)
+            y2=np.delete(y2,pos1[0],axis=0)
+            z1=np.delete(z1,pos1[0],axis=0)
+            z2=np.delete(z2,pos1[0],axis=0)
+            #next connection
+            pos1=np.where(edge_l==e1)
+            
+        #same thing the other way (to the right)
+        pos2=np.where(edge_l==e2)
+
+        while len(pos2[0])>0:
+            
+            if np.all(pos2[1]==0):
+                xc=np.vstack((xc,x2[pos2[0]]))
+                yc=np.vstack((yc,y2[pos2[0]]))
+                zc=np.vstack((zc,z2[pos2[0]]))
+                #next edge:
+                e2=edge_l[pos2[0],1]
+            else:
+                xc=np.vstack((xc,x1[pos2[0]]))
+                yc=np.vstack((yc,y1[pos2[0]]))
+                zc=np.vstack((zc,z1[pos2[0]]))
+                #next edge:
+                e2=edge_l[pos2[0],0]
+                
+            #erase the lines of this
+            edge_l=np.delete(edge_l,pos2[0],axis=0)
+            x1=np.delete(x1,pos2[0],axis=0)
+            x2=np.delete(x2,pos2[0],axis=0)
+            y1=np.delete(y1,pos2[0],axis=0)
+            y2=np.delete(y2,pos2[0],axis=0)
+            z1=np.delete(z1,pos2[0],axis=0)
+            z2=np.delete(z2,pos2[0],axis=0)
+            #next connection
+            pos2=np.where(edge_l==e2)
+            
+        #save xc,yc contour: 
+        newcontour = OrderedDict()
+        newcontour['nods'] = np.size(xc)
+        newcontour['density'] = 1 
+        newcontour['closed'] = 0
+        newcontour['x'] = np.ma.filled(xc.astype(float), np.nan)
+        newcontour['y'] = np.ma.filled(yc.astype(float), np.nan)
+        newcontour['z'] = np.ma.filled(zc.astype(float), np.nan)
+        newcontour['name'] = ''
+        contours.append(newcontour)
+        
+    return contours
Index: /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.py	(revision 26374)
+++ /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.py	(revision 26374)
@@ -0,0 +1,39 @@
+import numpy as np
+from model import model
+from contourlevelzero import *
+from ExpToLevelSet import *
+
+def reinitializelevelset(md,levelset):
+    """REINITIALIZELEVELSET - reinitialize levelset as a signed distance function
+
+    Usage:
+       levelsetnew = reinitializelevelset(md,levelset)
+    """
+    
+    # if md is 3d, levelset should be projected on a 2d mesh 
+    if len(levelset) == 0:
+        raise IOError("levelset provided is empty")
+    
+    if md.mesh.dimension() == 3:
+        if len(levelset)!=md.mesh.numberofvertices2d:
+            raise IOError("levelset provided should be specified at the 2d vertices of the mesh")
+        else:
+            if len(levelset)!=md.mesh.numberofvertices:
+                raise IOError("levelset provided should be specified at the vertices of the mesh")
+            
+    #First: extract segments
+    contours=contourlevelzero(md,levelset,0)
+    
+    #Now, make this a distance field (might not be closed)
+    levelsetnew=np.abs(ExpToLevelSet(md.mesh.x,md.mesh.y,contours)).T # levelsetnew comes on the 3d vertices, if mesh is 3d
+    
+    #Finally, change sign
+    pos=np.where(levelset<0) # if mesh is 3d, it refers to the vertices on the base
+    if md.mesh.dimension()==3:
+        for i in range(md.mesh.numberoflayers):
+            pos3d=pos[0]+i*md.mesh.numberofvertices2d
+            levelsetnew[pos3d]=-levelsetnew[pos3d]
+    else:
+        levelsetnew[pos[0]]=-1*levelsetnew[pos[0]]
+        
+    return levelsetnew