Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17502)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 17503)
@@ -4,4 +4,6 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../cores/cores.h"
 
 #include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
@@ -29,8 +31,12 @@
 		}
 	}
+
+	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 	int finiteelement=P1Enum;
+	if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 	::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -46,4 +52,11 @@
 	/*parameters: */
 	bool save_results;
+	int maxiter = 3;
+	int step;
+	IssmDouble reltol = 0.05;
+
+	Vector<IssmDouble>* lsfg     = NULL;
+	Vector<IssmDouble>* lsfg_old = NULL;
+
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
@@ -53,13 +66,36 @@
 	/* 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);
+	SetReinitSPCs(femmodel);
+
+	step = 1;
+	for(;;){
+
+		_printf_("smoothing lsf slope\n");
+		/* smoothen slope of lsf for computation of normal on ice domain*/
+		levelsetfunctionslope_core(femmodel);
+
+		//solve current artificial time step
+		if(VerboseSolution()) _printf0_("call computational core for reinitialization in step " << step << ":\n");
+		solutionsequence_linear(femmodel);
+		GetSolutionFromInputsx(&lsfg,femmodel);
+
+		if(step>1){
+			if(VerboseSolution()) _printf0_("   checking reinitialization convergence\n");
+			if(ReinitConvergence(lsfg,lsfg_old,reltol)) break;
+		}
+		if(step>maxiter){
+			if(VerboseSolution()) _printf0_("   maximum number reinitialization iterations " << maxiter << " reached\n");
+			break;
+		}
+
+		/*update results and increase counter*/
+		delete lsfg_old;lsfg_old=lsfg;
+		step++;
+	}
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		int outputs = MaskIceLevelsetEnum;
-		femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+		int outputs[1] = {MaskIceLevelsetEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
 	}
 
@@ -73,12 +109,11 @@
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
-
+	
 	/*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 Jdet,D_scalar;
+	IssmDouble dtau = 1.;
+	IssmDouble mu = 1.;
 	IssmDouble* xyz_list = NULL;
 
@@ -88,12 +123,11 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    basis    = xNew<IssmDouble>(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);
+	IssmDouble*    D		= xNew<IssmDouble>(dim*dim);
+	IssmDouble*    dlsf	= xNew<IssmDouble>(dim);
+	IssmDouble*    normal= xNew<IssmDouble>(dim);
+
 	element->GetVerticesCoordinates(&xyz_list);
-	h = element->CharacteristicLength();
 
 	/* Start  looping on the number of gaussian points: */
@@ -103,63 +137,39 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,xyz_list,gauss);
+		D_scalar=gauss->weight*Jdet;
+
+		if(dtau!=0.){
+			element->NodalFunctions(basis,gauss);
+			TripleMultiply(basis,numnodes,1,0,
+						&D_scalar,1,1,0,
+						basis,1,numnodes,0,
+						&Ke->values[0],1);
+			D_scalar*=dtau;
+		}
+
 		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];
+					D[row*dim+col]=mu*D_scalar;
 				else
-					D[row][col]=0.;
-		TripleMultiply(B,dim,numnodes,1,
-					&D[0][0],dim,dim,0,
+					D[row*dim+col]=0.;
+		TripleMultiply(Bprime,dim,numnodes,1,
+					D,dim,dim,0,
 					Bprime,dim,numnodes,0,
 					&Ke->values[0],1);
 
-		/* Stabilization *//*{{{*/
-		stabilization=1;
+		/* Stabilization */
+		stabilization=0;
 		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>(basis);
 	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
 	delete gauss;
 	return Ke;
@@ -168,15 +178,30 @@
 	
 	/*Intermediaries */
-	IssmDouble Jdet;
-	IssmDouble* xyz_list = NULL;
-	
+	int i,k;
+	int dim = 2;
+	IssmDouble dtau = 1.;
+	IssmDouble mu = 1.;
+	IssmDouble Jdet, D_scalar;
+	IssmDouble lsf;
+	IssmDouble norm_dlsf;
+	IssmDouble dbasis_normal;
+
 	/*Fetch number of nodes */
 	int numnodes = element->GetNumberOfNodes();
 
+	IssmDouble* xyz_list = NULL;
 	IssmDouble* basis = xNew<IssmDouble>(numnodes);
+	IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
+	IssmDouble* dlsf = xNew<IssmDouble>(dim);
+	IssmDouble* normal = xNew<IssmDouble>(dim);
 	element->GetVerticesCoordinates(&xyz_list);
 
 	/*Initialize Element vector*/
 	ElementVector* pe = element->NewElementVector();
+
+	/*Retrieve all inputs and parameters*/
+	Input* lsf_input = element->GetInput(MaskIceLevelsetEnum); _assert_(lsf_input);
+	Input* lsf_slopex_input=element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+	Input* lsf_slopey_input=element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
 
 	Gauss* gauss=element->NewGauss(2);
@@ -186,20 +211,101 @@
 		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]; 
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		D_scalar=Jdet*gauss->weight;
+
+		if(dtau!=0.){
+			/* old function value */
+			lsf_input->GetInputValue(&lsf,gauss);
+			for(i=0;i<numnodes;i++) pe->values[i]+=D_scalar*lsf*basis[i];
+			D_scalar*=dtau;
+		}
+
+		lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+		lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+
+		/*get normal*/
+		norm_dlsf=0.;
+		for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i]; 
+		norm_dlsf=sqrt(norm_dlsf);
+		if(norm_dlsf>0.)
+			for(i=0;i<dim;i++)	normal[i]=dlsf[i]/norm_dlsf;
+		else
+			for(i=0;i<dim;i++)	normal[i]=0.;
+
+		/* multiply normal and dbasis */
+		for(i=0;i<numnodes;i++){
+			dbasis_normal=0.;
+			for(k=0;k<dim;k++) dbasis_normal+=dbasis[k*numnodes+i]*normal[k];
+			pe->values[i]+=D_scalar*mu*dbasis_normal; 
+		}
 	}
 
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(dlsf);
+	xDelete<IssmDouble>(normal);
 	return pe;
 	}/*}}}*/
 void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	IssmDouble   lsf;
+	int          meshtype,dim,approximation,dofpernode;
+	int*         doflist = NULL;
+
+	/*Get some parameters*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum: dim = 2; dofpernode = 1; break;
+		case Mesh2DverticalEnum:   dim = 2; dofpernode = 1; break;
+		case Mesh3DEnum:           dim = 3; dofpernode = 1; break;
+		case Mesh3DtetrasEnum:     dim = 3; dofpernode = 1; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*dofpernode;
+
+	/*Fetch dof list and allocate solution vector*/
+	element->GetDofList(&doflist,approximation,GsetEnum);
+	IssmDouble* values = xNew<IssmDouble>(numdof);
+
+	/*Get inputs*/
+	Input* lsf_input=element->GetInput(MaskIceLevelsetEnum); _assert_(lsf_input);
+
+	Gauss* gauss=element->NewGauss();
+	for(int i=0;i<numnodes;i++){
+		gauss->GaussNode(element->FiniteElement(),i);
+
+		lsf_input->GetInputValue(&lsf,gauss);
+		values[i*dofpernode+0]=lsf;
+	}
+
+	solution->SetValues(numdof,doflist,values,INS_VAL);
+
+	/*Free ressources:*/
+	delete gauss;
+	xDelete<IssmDouble>(values);
+	xDelete<int>(doflist);
+
 }/*}}}*/
 void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	int meshtype;
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum);
+			break;
+		case Mesh3DEnum:
+			element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum);
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
 }/*}}}*/
 void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	_error_("not implemented yet");
+	/* Do nothing for now */
 }/*}}}*/
 void LsfReinitializationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
@@ -260,5 +366,5 @@
 
 /* Other */
-void LsfReinitializationAnalysis::UpdateReinitSPCs(FemModel* femmodel){/*{{{*/
+void LsfReinitializationAnalysis::SetReinitSPCs(FemModel* femmodel){/*{{{*/
 
 	int i,k, numnodes;
@@ -270,6 +376,8 @@
 		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		for(k=0;k<element->GetNumberOfNodes();k++){
-				node=element->GetNode(k);
+			node=element->GetNode(k);
+			if(node->IsActive()){
 				node->DofInFSet(0); 
+			}
 		}
 	}
@@ -287,5 +395,7 @@
 			for(k=0;k<numnodes;k++){
 				node=element->GetNode(k);
-				node->ApplyConstraint(1,lsf[k]);
+				if(node->IsActive()){
+					node->ApplyConstraint(1,lsf[k]);
+				}
 			}
 			xDelete<IssmDouble>(lsf);
@@ -298,23 +408,40 @@
 	/* Intermediaries */
 	int i,k;
+	IssmDouble dmaxp,dmaxm,val;
+	Element* element;
 
 	/*Initialize vector with number of vertices*/
 	int numvertices=femmodel->vertices->NumberOfVertices();
-	Element* element;
 
 	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))
-			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))
+		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);
 	}
 
@@ -327,4 +454,5 @@
 	delete vec_dist_zerolevelset;
 	delete dist_zerolevelset;
+
 }/*}}}*/
 void LsfReinitializationAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
@@ -334,13 +462,16 @@
 
 	/* Intermediaries */
-	const int dim=3;
+	int dim=3;
 	int i,d;
+	IssmDouble dist,lsf_old;
+
 	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* s0 = xNew<IssmDouble>(dim);
+	IssmDouble* s1 = xNew<IssmDouble>(dim);
+	IssmDouble* v = xNew<IssmDouble>(dim);
 	IssmDouble* xyz_list = NULL;
 	IssmDouble* xyz_list_zero = NULL;
@@ -355,7 +486,7 @@
 
 	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];
+	for(i=0;i<dim;i++){
+		s0[i]=xyz_list_zero[0+i];
+		s1[i]=xyz_list_zero[3+i];
 	}
 
@@ -363,5 +494,5 @@
 	for(i=0;i<numvertices;i++){
 		for(d=0;d<dim;d++)
-			v[d]=xyz_list[3*i+d];
+			v[d]=xyz_list[dim*i+d];
 		dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
 		signed_dist[i]=sign_lsf[i]*dist;
@@ -371,5 +502,6 @@
 	for(i=0;i<numvertices;i++){
 		vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-		if(xIsNan<IssmDouble>(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
+		/* 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);
 	}
@@ -378,4 +510,8 @@
 	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){/*{{{*/
@@ -403,3 +539,28 @@
 	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;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17502)
+++ /issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h	(revision 17503)
@@ -33,8 +33,9 @@
 	
 	/* Other */
-	void UpdateReinitSPCs(FemModel* femmodel);
+	void SetReinitSPCs(FemModel* femmodel);
 	void SetDistanceOnIntersectedElements(FemModel* femmodel);
 	void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);
 	IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
+	bool ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol);
 };
 #endif
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp	(revision 17502)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp	(revision 17503)
@@ -76,4 +76,5 @@
 		if(solution_enum==TransientSolutionEnum && analysis_enum==LevelsetAnalysisEnum && islevelset==false) continue;
 		if(solution_enum==TransientSolutionEnum && analysis_enum==ExtrapolationAnalysisEnum && islevelset==false) continue;
+		if(solution_enum==TransientSolutionEnum && analysis_enum==LsfReinitializationAnalysisEnum && islevelset==false) continue;