Index: /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.m	(revision 24317)
+++ /issm/trunk-jpl/src/m/parameterization/reinitializelevelset.m	(revision 24318)
@@ -5,4 +5,12 @@
 %      levelsetnew = reinitializelevelset(md,levelset)
 
+% if md is 3d, levelset should be projected on a 2d mesh 
+
+if isempty(levelset), error('levelset provided is empty'); end
+if dimension(md.mesh)==3,
+   if length(levelset)~=md.mesh.numberofvertices2d, error('levelset provided should be specified at the 2d vertices of the mesh'); end
+else
+   if length(levelset)~=md.mesh.numberofvertices, error('levelset provided should be specified at the vertices of the mesh'); end
+end
 
 %First: extract segments
@@ -10,7 +18,14 @@
 
 %Now, make this a distance field (might not be closed)
-levelsetnew=abs(ExpToLevelSet(md.mesh.x,md.mesh.y,contours));
+levelsetnew=abs(ExpToLevelSet(md.mesh.x,md.mesh.y,contours)); % levelsetnew comes on the 3d vertices, if mesh is 3d
 
 %Finally, change sign
-pos = find(levelset<0);
-levelsetnew(pos) = -levelsetnew(pos);
+pos=find(levelset<0); % if mesh is 3d, it refers to the vertices on the base
+if dimension(md.mesh)==3
+	for i=1:md.mesh.numberoflayers
+		pos3d=pos+(i-1)*md.mesh.numberofvertices2d;
+		levelsetnew(pos3d)=-levelsetnew(pos3d);
+	end	
+else
+	levelsetnew(pos)=-levelsetnew(pos);
+end