Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17306)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17307)
@@ -6,47 +6,156 @@
 
 /*Model processing*/
-int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){
-	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){
-	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){
-	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){
-	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){
-	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){
-	_error_("not implemented yet");
-}
+int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
+	return 1;
+}/*}}}*/
+void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
 
 /*Finite element Analysis*/
-void           LsfReinitializationAnalysis::Core(FemModel* femmodel){
+void           LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
 	_error_("not implemented yet");
-}
-ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){
+}/*}}}*/
+ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){
+}/*}}}*/
+ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){
+}/*}}}*/
+ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){
+}/*}}}*/
+ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){
+}/*}}}*/
+void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){
+}/*}}}*/
+void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){
+}/*}}}*/
+void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	_error_("not implemented yet");
-}
+}/*}}}*/
 
+/* Other */
+void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+
+	/* Intermediaries */
+	int i;
+
+	/*Initialize vector with number of vertices*/
+	int numvertices=femmodel->vertices->NumberOfVertices();
+	Vector<IssmDouble>* vec_dist_zerolevelset=new Vector<IssmDouble>(numvertices); //vertices that have ice at next time step
+
+	/*TODO: Fill vector with values of old level set function: */
+	for(i=0;i<numvertices;i++){
+		vec_dist_zerolevelset->SetValue(i,0.,INS_VAL);
+	}
+	for(i=0;i<femmodel->elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(element->IsZeroLevelset(MaskIceLevelsetEnum))
+			SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
+	}
+
+	/*Assemble vector and serialize */
+	vec_dist_zerolevelset->Assemble();
+	IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
+	InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
+
+	/*Clean up and return*/
+	delete vec_dist_zerolevelset;
+	delete dist_zerolevelset;
+}/*}}}*/
+void LsfReinitializationAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
+
+	if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
+		return;
+
+	/* Intermediaries */
+	const int dim=3;
+	int i,d;
+	int numvertices=element->GetNumberOfVertices();
+	IssmDouble s0[dim], s1[dim], v[dim];
+	IssmDouble dist,lsf_old;
+
+	IssmDouble* lsf = xNew<IssmDouble>(numvertices);
+	IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
+	IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
+	IssmDouble* xyz_list = NULL;
+	IssmDouble* xyz_list_zero = NULL;
+
+	/* retrieve inputs and parameters */
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum);
+
+	/* get sign of levelset function */
+	for(i=0;i<numvertices;i++){
+		if(lsf[i]==0.)
+			sign_lsf[i]=1.;
+		else
+			sign_lsf[i]=lsf[i]/fabs(lsf[i]);
+	}
+
+	element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
+	for(d=0;d<dim;d++){
+		s0[d]=xyz_list_zero[0+d];
+		s1[d]=xyz_list_zero[3+d];
+	}
+
+	/* get signed_distance of vertices to zero levelset straight */
+	for(i=0;i<numvertices;i++){
+		for(d=0;d<dim;d++)
+			v[d]=xyz_list[3*i+d];
+		dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
+		signed_dist[i]=sign_lsf[i]*dist;
+	}
+	
+	/* insert signed_distance into vec_signed_dist, if computed distance is smaller */
+	for(i=0;i<numvertices;i++){
+		vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
+		if(fabs(signed_dist[i])<fabs(lsf_old))
+			vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
+	}
+
+	xDelete<IssmDouble>(lsf);
+	xDelete<IssmDouble>(sign_lsf);
+	xDelete<IssmDouble>(signed_dist);
+}/*}}}*/
+IssmDouble LsfReinitializationAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
+	// returns distance d of point q to straight going through points s0, s1
+	// d=|a x b|/|b|
+	// with a=q-s0, b=s1-s0
+	
+	/* Intermediaries */
+	const int dim=2;
+	int i;
+	IssmDouble a[dim], b[dim];
+	IssmDouble norm_b;
+
+	for(i=0;i<dim;i++){
+		a[i]=q[i]-s0[i];
+		b[i]=s1[i]-s0[i];
+	}
+	
+	norm_b=0.;
+	for(i=0;i<dim;i++)
+		norm_b+=b[i]*b[i];
+	norm_b=sqrt(norm_b);
+	_assert_(norm_b>0.);
+
+	return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
+}/*}}}*/
+
Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17306)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17307)
@@ -29,4 +29,9 @@
 	void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 	void UpdateConstraints(FemModel* femmodel);
+	
+	/* Other */
+	void SetDistanceOnIntersectedElements(FemModel* femmodel);
+	void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);
+	IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
 };
 #endif