Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19475)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19476)
@@ -143,4 +143,5 @@
 	basalelement->FindParam(&calvinglaw,CalvingLawEnum);
 	basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
+	calvinglaw=DefaultCalvingEnum;
 	switch(domaintype){
 		case Domain2DverticalEnum:   dim = 1; break;
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 19475)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 19476)
@@ -178,4 +178,6 @@
 		virtual void       CalvingRatePi(void)=0;
 		virtual void       CalvingRateDev(void){_error_("not implemented yet");};
+		virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
+		virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};
 		virtual IssmDouble CharacteristicLength(void)=0;
 		virtual void       ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19475)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19476)
@@ -350,5 +350,5 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 	IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-	IssmDouble  sigma_vm,sigma_max;
+	IssmDouble  sigma_vm,sigma_max,epse_2;
 
 	/* Get node coordinates and dof list: */
@@ -384,6 +384,13 @@
 
 		/*Calculate sigma_vm*/
-		sigma_vm  = sqrt(3./2.) * B * pow(lambda1*lambda1 + lambda2*lambda2,1./(2.*n));
+		epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
+		sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+		//sigma_max = 125.e+3;
 		sigma_max = 350.e+3;
+		sigma_max = 450.e+3;
+		sigma_max = 800.e+3; //too much
+		sigma_max = 700.e+3;
+		sigma_max = 670.e+3;
+		//sigma_max = 550.e+3;
 
 		/*Assign values*/
@@ -400,4 +407,124 @@
 	/*Clean up and return*/
 	delete gauss;
+}
+/*}}}*/
+void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
+
+	if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;
+
+	IssmDouble* xyz_list_zero = NULL;
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+	this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0], MaskIceLevelsetEnum);
+	if(xyz_list_zero){
+		IssmDouble x[2];
+		IssmDouble y[2];
+		x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0];
+		y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1];
+		segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
+	}
+	xDelete<IssmDouble>(xyz_list_zero);
+
+//	IssmDouble ls[NUMVERTICES];
+//	IssmDouble  xyz_list[NUMVERTICES][3];
+//
+//	if(IsIceInElement()){
+//
+//		/*Retrieve all inputs and parameters*/
+//		GetInputListOnVertices(&ls[0],levelset_enum);
+//
+//		/*If the level set is awlays <0, there is no ice front here*/
+//		bool iszerols= false;
+//		if(IsIceInElement()){
+//			if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
+//				iszerols = true;
+//			}
+//		}
+//
+//		if(iszerols){
+//			/*OK we have one segment!*/
+//			::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+//			int count = 0;
+//			IssmDouble x[2];
+//			IssmDouble y[2];
+//
+//			for(int i=0;i<NUMVERTICES,i++){
+//				int index1 = i;
+//				int index1 = (i+1)%3;
+//				if(ls[index1]<=0 && ls[index2]>=0){
+//
+//				}
+//
+//			}
+//			Contour* segment = new Contour<IssmDouble>(segment->Size()+1,2,x,y,false);
+//		}
+//
+//	}
+//
+//	_error_("STOP");
+}
+/*}}}*/
+void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/
+
+	/*Intermediaries*/
+	double 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);
+	InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
+
+	/*Get distance from list of segments and reset ls*/
+	for(int j=0;j<NUMVERTICES;j++){
+		double 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*/
+			double 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] = dmin;
+		}
+		else{
+			ls[j] = - dmin;
+		}
+	}
+
+	/*Update Levelset*/
+	this->inputs->AddInput(new TriaInput(MaskIceLevelsetEnum,&ls[0],P1Enum));
 }
 /*}}}*/
@@ -3144,5 +3271,5 @@
 	IssmDouble* xyz_zero=NULL;
 
-	this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+	this->GetLevelsetIntersection(&indices, &numiceverts, s,levelsetenum,0.);
 	
 	//TODO: check if for 2 iceverts front segment is oriented in CCW way
@@ -3165,5 +3292,5 @@
 	else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
 		IssmDouble lsf[NUMVERTICES];
-		this->GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+		this->GetInputListOnVertices(&lsf[0],levelsetenum);
 		counter=0;
 		for(i=0;i<NUMVERTICES;i++){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 19475)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 19476)
@@ -55,4 +55,6 @@
 		void			CalvingRatePi();
 		void			CalvingRateDev();
+		void			WriteLevelsetSegment(DataSet* segments);
+		void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
 		IssmDouble  CharacteristicLength(void);
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19475)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19476)
@@ -842,4 +842,55 @@
 		(element->*function)();
 	}
+
+}
+/*}}}*/
+void FemModel::ResetLevelset(void){/*{{{*/
+
+	/*recover my_rank:*/
+	int my_rank   = IssmComm::GetRank();
+	int num_procs = IssmComm::GetSize();
+
+	/*1: go throug all elements of this partition and figure out how many
+	 * segments we have (corresopnding to levelset = 0)*/
+	DataSet* segments=new DataSet();
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->WriteLevelsetSegment(segments);
+	}
+
+	/*2: now get the segments from all partitions*/
+	int  segcount=segments->Size();
+	int* allsegcount=xNew<int>(num_procs);
+	ISSM_MPI_Gather(&segcount,1,ISSM_MPI_INT,allsegcount,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(allsegcount,num_procs,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	/* Every cpu should start its own dof count at the end of the dofcount from cpu-1*/
+	int numseg_offset=0;
+	int numseg=0;
+	for(int i=0;i<my_rank;  i++) numseg_offset+=allsegcount[i];
+	for(int i=0;i<num_procs;i++) numseg+=allsegcount[i];
+	IssmDouble* segmentlist    = xNewZeroInit<IssmDouble>(4*numseg);
+	IssmDouble* allsegmentlist = xNewZeroInit<IssmDouble>(4*numseg);
+	for(int i=0;i<segments->Size();i++){
+		Contour<IssmDouble>* segment=(Contour<IssmDouble>*)segments->GetObjectByOffset(i);
+		_assert_(segment->nods == 2);
+		segmentlist[(numseg_offset+i)*4 + 0] = segment->x[0];
+		segmentlist[(numseg_offset+i)*4 + 1] = segment->y[0];
+		segmentlist[(numseg_offset+i)*4 + 2] = segment->x[1];
+		segmentlist[(numseg_offset+i)*4 + 3] = segment->y[1];
+	}
+	ISSM_MPI_Allreduce((void*)segmentlist,(void*)allsegmentlist,4*numseg,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+	delete segments;
+	xDelete<IssmDouble>(segmentlist);
+	xDelete<int>(allsegcount);
+
+	/*3: update levelset for all elements*/
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->ResetLevelsetFromSegmentlist(allsegmentlist,numseg);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(allsegmentlist);
 
 }
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 19475)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 19476)
@@ -94,4 +94,5 @@
 		void CalvingRatePix();
 		void CalvingRateDevx();
+		void ResetLevelset();
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);