Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17333)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17334)
@@ -5,4 +5,6 @@
 #include "../modules/modules.h"
 
+#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
+
 /*Model processing*/
 int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
@@ -13,8 +15,22 @@
 }/*}}}*/
 void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-	/* Do nothing for now */
+	int    finiteelement;
+
+	/*Finite element type*/
+	finiteelement = P1Enum;
+
+	/*Update elements: */
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
+			counter++;
+		}
+	}
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	/* Do nothing for now */
+	int finiteelement=P1Enum;
+	::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -26,9 +42,30 @@
 
 /*Finite element Analysis*/
-void           LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
-	_error_("not implemented yet");
+void  LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
+
+	/*parameters: */
+	bool save_results;
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	/*activate formulation: */
+	femmodel->SetCurrentConfiguration(LsfReinitializationAnalysisEnum);
+
+	/* set spcs for reinitialization */
+	if(VerboseSolution()) _printf0_("Update spcs for reinitialization:\n");
+	UpdateReinitSPCs(femmodel);
+
+	if(VerboseSolution()) _printf0_("call computational core for reinitialization:\n");
+// 	solutionsequence_lsfreinit_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs = MaskIceLevelsetEnum;
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+	}
+
 }/*}}}*/
 ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){/*{{{*/
-	_error_("not implemented yet");
+	/*Default, return NULL*/
+	return NULL;
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
@@ -36,9 +73,125 @@
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries */
+	const int dim = 2;
+	int        i,row,col,stabilization;
+	IssmDouble Jdet,D_scalar,h;
+	IssmDouble dlsf[3],normal[3];
+	IssmDouble norm_dlsf;
+	IssmDouble hx,hy,hz,kappa;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble     D[dim][dim];
+
+	/*Retrieve all inputs and parameters*/
+	Input* lsfpicard_input=element->GetInput(LevelsetfunctionPicardEnum); _assert_(lsfpicard_input);
+	element->GetVerticesCoordinates(&xyz_list);
+	h = element->CharacteristicLength();
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		GetB(B,element,xyz_list,gauss);
+		GetBprime(Bprime,element,xyz_list,gauss);
+
+		/* Get normal from last iteration on lsf */
+		lsfpicard_input->GetInputDerivativeValue(&dlsf[0],xyz_list,gauss);
+
+		norm_dlsf=0.;
+		for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i]; 
+		norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
+		for(i=0;i<dim;i++)
+			normal[i]=dlsf[i]/norm_dlsf;
+		
+		D_scalar=gauss->weight*Jdet;
+
+		for(row=0;row<dim;row++)
+			for(col=0;col<dim;col++)
+				if(row==col)
+					D[row][col]=D_scalar*normal[row];
+				else
+					D[row][col]=0.;
+		TripleMultiply(B,dim,numnodes,1,
+					&D[0][0],dim,dim,0,
+					Bprime,dim,numnodes,0,
+					&Ke->values[0],1);
+
+		/* Stabilization *//*{{{*/
+		stabilization=1;
+		if (stabilization==0){/* no stabilization, do nothing*/}
+		else if(stabilization==1){
+			/* Artificial Diffusion */
+			element->ElementSizes(&hx,&hy,&hz);
+			h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+			kappa=h/2.; 
+			D[0][0]=D_scalar*kappa;
+			D[0][1]=0.;
+			D[1][0]=0.;
+			D[1][1]=D_scalar*kappa;
+			TripleMultiply(Bprime,dim,numnodes,1,
+						&D[0][0],dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		}
+		else if(stabilization==2){
+			/*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
+			for(row=0;row<dim;row++)
+				for(col=0;col<dim;col++)
+					D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+
+			TripleMultiply(Bprime,dim,numnodes,1,
+						&D[0][0],dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		}/*}}}*/
+	}/*}}}*/
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
 }/*}}}*/
 ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
+	
+	/*Intermediaries */
+	IssmDouble Jdet;
+	IssmDouble* xyz_list = NULL;
+	
+	/*Fetch number of nodes */
+	int numnodes = element->GetNumberOfNodes();
+
+	IssmDouble* basis = xNew<IssmDouble>(numnodes);
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/*Initialize Element vector*/
+	ElementVector* pe = element->NewElementVector();
+
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*basis[i]; 
+	}
+
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	return pe;
+	}/*}}}*/
 void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
@@ -50,21 +203,116 @@
 	_error_("not implemented yet");
 }/*}}}*/
+void LsfReinitializationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ N ]
+	 *          [ N ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = basis[i];
+		B[numnodes*1+i] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void LsfReinitializationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ dN/dx ]
+	 *                [ dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+
+}/*}}}*/
 
 /* Other */
+void LsfReinitializationAnalysis::UpdateReinitSPCs(FemModel* femmodel){/*{{{*/
+
+	int i,k, numnodes;
+	Element* element;
+	Node* node;
+
+	/* deactivate all spcs */
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		for(k=0;k<element->GetNumberOfNodes();k++){
+				node=element->GetNode(k);
+				node->DofInFSet(0); 
+		}
+	}
+
+	SetDistanceOnIntersectedElements(femmodel);
+
+	/* reactivate spcs on elements intersected by zero levelset */
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
+			/*iterate over nodes and set spc */
+			numnodes=element->GetNumberOfNodes();
+			IssmDouble* lsf = xNew<IssmDouble>(numnodes);
+			element->GetInputListOnNodes(&lsf[0],MaskIceLevelsetEnum);
+			for(k=0;k<numnodes;k++){
+				node=element->GetNode(k);
+				node->ApplyConstraint(1,lsf[k]);
+			}
+			xDelete<IssmDouble>(lsf);
+		}
+	}
+
+}/*}}}*/
 void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
 
 	/* Intermediaries */
-	int i;
+	int i,k;
 
 	/*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);
-	}
+	Element* element;
+
+	Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
+	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
+	
 	for(i=0;i<femmodel->elements->Size();i++){
-		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(element->IsZeroLevelset(MaskIceLevelsetEnum))
+			for(k=0;k<element->GetNumberOfVertices();k++)
+				vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL); 
+	}
+
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		if(element->IsZeroLevelset(MaskIceLevelsetEnum))
 			SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
@@ -103,10 +351,6 @@
 
 	/* 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]);
-	}
+	for(i=0;i<numvertices;i++)
+		sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
 
 	element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
@@ -127,5 +371,5 @@
 	for(i=0;i<numvertices;i++){
 		vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-		if(fabs(signed_dist[i])<fabs(lsf_old))
+		if(isnan(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
 			vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
 	}
Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17333)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17334)
@@ -29,6 +29,9 @@
 	void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 	void UpdateConstraints(FemModel* femmodel);
+	void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+	void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 	
 	/* Other */
+	void UpdateReinitSPCs(FemModel* femmodel);
 	void SetDistanceOnIntersectedElements(FemModel* femmodel);
 	void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17333)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17334)
@@ -434,4 +434,5 @@
 				name==LevelsetfunctionSlopeXEnum ||
 				name==LevelsetfunctionSlopeYEnum ||
+				name==LevelsetfunctionPicardEnum ||
 				name==GradientEnum ||
 				name==OldGradientEnum  ||
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17333)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17334)
@@ -682,4 +682,5 @@
 	LevelsetfunctionSlopeXEnum,
 	LevelsetfunctionSlopeYEnum,
+	LevelsetfunctionPicardEnum,
 	/*}}}*/
 	MaximumNumberOfDefinitionsEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17333)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17334)
@@ -641,4 +641,5 @@
 		case LevelsetfunctionSlopeXEnum : return "LevelsetfunctionSlopeX";
 		case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY";
+		case LevelsetfunctionPicardEnum : return "LevelsetfunctionPicard";
 		case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
 		default : return "unknown";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17333)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17334)
@@ -656,4 +656,5 @@
 	      else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
+	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
 	      else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=7;
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17333)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17334)
@@ -633,3 +633,4 @@
 def LevelsetfunctionSlopeXEnum(): return StringToEnum("LevelsetfunctionSlopeX")[0]
 def LevelsetfunctionSlopeYEnum(): return StringToEnum("LevelsetfunctionSlopeY")[0]
+def LevelsetfunctionPicardEnum(): return StringToEnum("LevelsetfunctionPicard")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]
Index: /issm/trunk-jpl/src/m/enum/LevelsetfunctionPicardEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/LevelsetfunctionPicardEnum.m	(revision 17334)
+++ /issm/trunk-jpl/src/m/enum/LevelsetfunctionPicardEnum.m	(revision 17334)
@@ -0,0 +1,11 @@
+function macro=LevelsetfunctionPicardEnum()
+%LEVELSETFUNCTIONPICARDENUM - Enum of LevelsetfunctionPicard
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=LevelsetfunctionPicardEnum()
+
+macro=StringToEnum('LevelsetfunctionPicard');