Index: issm/trunk/src/c/analyses/LsfReinitializationAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 18301)
+++ issm/trunk/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 19105)
@@ -10,9 +10,18 @@
 
 /*Model processing*/
+void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+	int finiteelement=P1Enum;
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+	::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+}/*}}}*/
 int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int domaintype,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){/*{{{*/
@@ -34,19 +43,10 @@
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 }/*}}}*/
-void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	int finiteelement=P1Enum;
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-	::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
-	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}/*}}}*/
-void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 	/* 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){/*{{{*/
 
 	/*parameters: */
@@ -249,5 +249,84 @@
 	return pe;
 	}/*}}}*/
-void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+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);
+
+}/*}}}*/
+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;
+}/*}}}*/
+void           LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
 	IssmDouble   lsf;
@@ -284,8 +363,8 @@
 
 }/*}}}*/
-void LsfReinitializationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           LsfReinitializationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int domaintype;
@@ -301,65 +380,142 @@
 	}
 }/*}}}*/
-void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	/* Do nothing for now */
-}/*}}}*/
-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::SetReinitSPCs(FemModel* femmodel){/*{{{*/
+bool           LsfReinitializationAnalysis::ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol){/*{{{*/
+
+	/*Output*/
+	bool converged = true;
+
+	/*Intermediary*/
+	Vector<IssmDouble>* dlsfg    = NULL;
+	IssmDouble          norm_dlsf,norm_lsf;
+
+	/*compute norm(du)/norm(u)*/
+	dlsfg=lsfg_old->Duplicate(); lsfg_old->Copy(dlsfg); dlsfg->AYPX(lsfg,-1.0);
+	norm_dlsf=dlsfg->Norm(NORM_TWO); norm_lsf=lsfg_old->Norm(NORM_TWO);
+	if (xIsNan<IssmDouble>(norm_dlsf) || xIsNan<IssmDouble>(norm_lsf)) _error_("convergence criterion is NaN!");
+	if((norm_dlsf/norm_lsf)<reltol){
+		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" < "<<reltol*100<<" %\n");
+	}
+	else{ 
+		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" > "<<reltol*100<<" %\n");
+		converged=false;
+	}
+
+	/*Cleanup*/
+	delete dlsfg;
+
+	return converged;
+}/*}}}*/
+void           LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+
+	/* Intermediaries */
+	int i;
+	IssmDouble dmaxp,dmaxm,val;
+	Element* element;
+
+	/*Initialize vector with number of vertices*/
+	int numvertices=femmodel->vertices->NumberOfVertices();
+
+	Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
+	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
+	
+	/* set distance on elements intersected by zero levelset */
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
+			SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
+		}
+	}
+	vec_dist_zerolevelset->Assemble();
+
+	/* Get maximum distance to interface along vertices */
+	dmaxp=0.; dmaxm=0.;
+	for(i=0;i<numvertices;i++){
+		vec_dist_zerolevelset->GetValue(&val,i); 
+		if((val>0.) && (val>dmaxp))
+			 dmaxp=val;
+		else if((val<0.) && (val<dmaxm))
+			 dmaxm=val;
+	}
+	//wait until all values are computed
+
+	/* set all none intersected vertices to max/min distance */
+	for(i=0;i<numvertices;i++){
+		vec_dist_zerolevelset->GetValue(&val,i);
+		if(val==1.) //FIXME: improve check
+			vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
+		else if(val==-1.)
+			vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
+	}
+
+	/*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 */
+	int dim=3;
+	int i,d;
+	IssmDouble dist,lsf_old;
+
+	int numvertices=element->GetNumberOfVertices();
+
+	IssmDouble* lsf = xNew<IssmDouble>(numvertices);
+	IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
+	IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
+	IssmDouble* s0 = xNew<IssmDouble>(dim);
+	IssmDouble* s1 = xNew<IssmDouble>(dim);
+	IssmDouble* v = xNew<IssmDouble>(dim);
+	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++)
+		sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
+
+	element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
+	for(i=0;i<dim;i++){
+		s0[i]=xyz_list_zero[0+i];
+		s1[i]=xyz_list_zero[3+i];
+	}
+
+	/* 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[dim*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());
+		/* initial lsf values are +-1. Update those fields or if distance to interface smaller.*/
+		if(fabs(lsf_old)==1. || 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);
+	xDelete<IssmDouble>(s0);
+	xDelete<IssmDouble>(s1);
+	xDelete<IssmDouble>(v);
+
+}/*}}}*/
+void           LsfReinitializationAnalysis::SetReinitSPCs(FemModel* femmodel){/*{{{*/
 
 	int i,k, numnodes;
@@ -369,5 +525,5 @@
 	/* deactivate all spcs */
 	for(i=0;i<femmodel->elements->Size();i++){
-		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		for(k=0;k<element->GetNumberOfNodes();k++){
 			node=element->GetNode(k);
@@ -382,5 +538,5 @@
 	/* reactivate spcs on elements intersected by zero levelset */
 	for(i=0;i<femmodel->elements->Size();i++){
-		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
 			/*iterate over nodes and set spc */
@@ -399,163 +555,5 @@
 
 }/*}}}*/
-void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
-
-	/* Intermediaries */
-	int i;
-	IssmDouble dmaxp,dmaxm,val;
-	Element* element;
-
-	/*Initialize vector with number of vertices*/
-	int numvertices=femmodel->vertices->NumberOfVertices();
-
-	Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
-	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
-	
-	/* set distance 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)){
-			SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
-		}
-	}
-	vec_dist_zerolevelset->Assemble();
-
-	/* Get maximum distance to interface along vertices */
-	dmaxp=0.; dmaxm=0.;
-	for(i=0;i<numvertices;i++){
-		vec_dist_zerolevelset->GetValue(&val,i); 
-		if((val>0.) && (val>dmaxp))
-			 dmaxp=val;
-		else if((val<0.) && (val<dmaxm))
-			 dmaxm=val;
-	}
-	//wait until all values are computed
-
-	/* set all none intersected vertices to max/min distance */
-	for(i=0;i<numvertices;i++){
-		vec_dist_zerolevelset->GetValue(&val,i);
-		if(val==1.) //FIXME: improve check
-			vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
-		else if(val==-1.)
-			vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
-	}
-
-	/*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 */
-	int dim=3;
-	int i,d;
-	IssmDouble dist,lsf_old;
-
-	int numvertices=element->GetNumberOfVertices();
-
-	IssmDouble* lsf = xNew<IssmDouble>(numvertices);
-	IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
-	IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
-	IssmDouble* s0 = xNew<IssmDouble>(dim);
-	IssmDouble* s1 = xNew<IssmDouble>(dim);
-	IssmDouble* v = xNew<IssmDouble>(dim);
-	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++)
-		sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
-
-	element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
-	for(i=0;i<dim;i++){
-		s0[i]=xyz_list_zero[0+i];
-		s1[i]=xyz_list_zero[3+i];
-	}
-
-	/* 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[dim*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());
-		/* initial lsf values are +-1. Update those fields or if distance to interface smaller.*/
-		if(fabs(lsf_old)==1. || 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);
-	xDelete<IssmDouble>(s0);
-	xDelete<IssmDouble>(s1);
-	xDelete<IssmDouble>(v);
-
-}/*}}}*/
-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;
-}/*}}}*/
-bool LsfReinitializationAnalysis::ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol){/*{{{*/
-
-	/*Output*/
-	bool converged = true;
-
-	/*Intermediary*/
-	Vector<IssmDouble>* dlsfg    = NULL;
-	IssmDouble          norm_dlsf,norm_lsf;
-
-	/*compute norm(du)/norm(u)*/
-	dlsfg=lsfg_old->Duplicate(); lsfg_old->Copy(dlsfg); dlsfg->AYPX(lsfg,-1.0);
-	norm_dlsf=dlsfg->Norm(NORM_TWO); norm_lsf=lsfg_old->Norm(NORM_TWO);
-	if (xIsNan<IssmDouble>(norm_dlsf) || xIsNan<IssmDouble>(norm_lsf)) _error_("convergence criterion is NaN!");
-	if((norm_dlsf/norm_lsf)<reltol){
-		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" < "<<reltol*100<<" %\n");
-	}
-	else{ 
-		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" > "<<reltol*100<<" %\n");
-		converged=false;
-	}
-
-	/*Cleanup*/
-	delete dlsfg;
-
-	return converged;
-}/*}}}*/
+void           LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/* Do nothing for now */
+}/*}}}*/
