Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22987)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22988)
@@ -3301,19 +3301,60 @@
 }
 /*}}}*/
-void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* distances){/*{{{*/
-
+void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble d,xn,yn;
 
 	/*Get current levelset and vertex coordinates*/
 	IssmDouble ls[NUMVERTICES];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
 	/*Get distance from list of segments and reset ls*/
 	for(int j=0;j<NUMVERTICES;j++){
+		IssmDouble dmin = 1.e+50;
+		for(int i=0;i<numsegments;i++){
+			IssmDouble x = xyz_list[j][0];
+			IssmDouble y = xyz_list[j][1];
+			IssmDouble l2 = (segments[4*i+2]-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (segments[4*i+3]-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]);
+
+			/*Segment has a length of 0*/
+			if(l2==0.){
+				d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+				if(d<dmin) dmin = d;
+				continue;
+			}
+
+			/*Consider the line extending the segment, parameterized as v + t (w - v).
+			 *We find projection of point p onto the line.
+			 *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+			IssmDouble t = ((x-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (y-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]))/l2;
+			if(t < 0.0){
+				// Beyond the 'v' end of the segment
+				d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+			}
+			else if (t > 1.0){
+				// Beyond the 'w' end of the segment
+				d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]);
+			}
+			else{
+				// Projection falls on the segment
+				xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]);
+				yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]);
+				d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
+			}
+
+			if(d<dmin) dmin = d;
+		}
+
+		/*Update signed distance*/
+		dmin = sqrt(dmin);
 		if(dmin>10000) dmin=10000;
 		if(ls[j]>0){
-			ls[j] = distances[this->vertices->Pid()];
+			ls[j] = dmin;
 		}
 		else{
-			ls[j] = - distances[this->vertices->Pid()];
+			ls[j] = - dmin;
 		}
 	}
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22987)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22988)
@@ -1050,54 +1050,24 @@
 
 	/*3: Add distance input to all elements*/
-	IssmDouble* distances = xNew<IssmDouble>(vertices->Size());
-	IssmDouble  d,xn,yn;
-	for(int i=0;i<vertices->Size();i++){
-		Vertex* vertex=dynamic_cast<Vertex*>(this->vertices->GetObjectByOffset(i));
-		IssmDouble dmin = 1.e+50;
-
-		for(int i=0;i<numseg;i++){
-			IssmDouble x = vertex->x;
-			IssmDouble y = vertex->y;
-			IssmDouble l2 = (allsegmentlist[4*i+2]-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (allsegmentlist[4*i+3]-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]);
-
-			/*Segment has a length of 0*/
-			if(l2==0.){
-				d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
-				if(d<dmin) dmin = d;
-				continue;
-			}
-
-			/*Consider the line extending the segment, parameterized as v + t (w - v).
-			 *We find projection of point p onto the line.
-			 *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
-			IssmDouble t = ((x-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (y-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]))/l2;
-			if(t < 0.0){
-				// Beyond the 'v' end of the segment
-				d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
-			}
-			else if (t > 1.0){
-				// Beyond the 'w' end of the segment
-				d = (x-allsegmentlist[4*i+2])*(x-allsegmentlist[4*i+2])+(y-allsegmentlist[4*i+3])*(y-allsegmentlist[4*i+3]);
-			}
-			else{
-				// Projection falls on the segment
-				xn = allsegmentlist[4*i+0] + t * (allsegmentlist[4*i+2] - allsegmentlist[4*i+0]);
-				yn = allsegmentlist[4*i+1] + t * (allsegmentlist[4*i+3] - allsegmentlist[4*i+1]);
-				d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
-			}
-
-			if(d<dmin) dmin = d;
-		}
-
-		/*Update signed distance*/
-		distances[vertex->pid] = sqrt(dmin);
-	}
-
-	//element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum);
-	InputUpdateFromVectorx(this,distances,distanceenum,VertexPIdEnum);
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		if(!element->IsOnBase()) continue;
+		element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum);
+	}
+
+	/*Extrude if necessary*/
+	int elementtype;
+	this->parameters->FindParam(&elementtype,MeshElementtypeEnum);
+	if(elementtype==PentaEnum){
+		InputExtrudex(this,distanceenum,-1);
+	}
+	else if(elementtype==TriaEnum){
+		/*no need to extrude*/
+	}
+	else{
+		_error_("not implemented yet");
+	}
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(distances);
 	xDelete<IssmDouble>(allsegmentlist);
 }/*}}}*/