Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18354)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18355)
@@ -194,11 +194,9 @@
 
 }/*}}}*/
-void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble deltat){/*{{{*/
+void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble* ulmin,IssmDouble* ulmax,IssmDouble deltat){/*{{{*/
 
 	/*Intermediaries*/
 	Vec         Ri_plus  = NULL;
 	Vec         Ri_minus = NULL;
-	IssmDouble  uiLmax   = 3.;
-	IssmDouble  uiLmin   = 2.;
 	int         ncols,ncols2,rstart,rend;
 	double      d;
@@ -234,6 +232,8 @@
 
 		/*Compute Qis and Ris*/
-		double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u[row]);
-		double Qi_minus = ml_serial[row]/deltat*(uiLmin - u[row]);
+		double Qi_plus  = ml_serial[row]/deltat*(3. - u[row]);
+		double Qi_minus = ml_serial[row]/deltat*(2. - u[row]);
+		//double Qi_plus  = ml_serial[row]/deltat*(ulmax[row] - u[row]);
+		//double Qi_minus = ml_serial[row]/deltat*(ulmin[row] - u[row]);
 		d = 1.;
 		if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
@@ -303,118 +303,123 @@
 	*pFbar = Fbar;
 }/*}}}*/
+void UpdateSolution(Vec u,Vec udot,Vec Ml,Vec Fbar,IssmDouble deltat){/*{{{*/
+
+	/*Intermediary*/
+	Vec temp1 = NULL;
+
+	/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
+	VecDuplicate(u,&temp1);
+	VecPointwiseDivide(temp1,Fbar,Ml); //temp1 = Ml^-1 temp2
+	VecAXPBY(udot,1.,1.,temp1);
+	VecAXPY(u,deltat,temp1);
+
+	/*CLean up and return*/
+	VecFree(&temp1);
+}/*}}}*/
 #endif
-	void solutionsequence_fct(FemModel* femmodel){
-
-		/*intermediary: */
-		Vector<IssmDouble>*  Ml = NULL;
-		Matrix<IssmDouble>*  K  = NULL;
-		Matrix<IssmDouble>*  Mc = NULL;
-		Vector<IssmDouble>*  ug = NULL;
-		Vector<IssmDouble>*  uf = NULL;
-
-		IssmDouble theta,deltat,dmax;
-		int        dof,ncols,ncols2,rstart,rend;
-		int        configuration_type;
-		double     d;
-		int*       cols  = NULL;
-		int*       cols2 = NULL;
-		double*    vals  = NULL;
-		double*    vals2 = NULL;
-
-		/*Create analysis*/
-		MasstransportAnalysis* analysis = new MasstransportAnalysis();
-
-		/*Recover parameters: */
-		femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
-		femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-		femmodel->UpdateConstraintsx();
-		theta = 0.5;
-
-		/*Create lumped mass matrix*/
-		analysis->LumpedMassMatrix(&Ml,femmodel);
-		analysis->MassMatrix(&Mc,femmodel);
-		analysis->FctKMatrix(&K,NULL,femmodel);
-
-		/*Convert matrices to PETSc matrices*/
-		Mat D_petsc  = NULL;
-		Mat LHS      = NULL;
-		Vec RHS      = NULL;
-		Vec u        = NULL;
-		Vec udot     = NULL;
-		Mat K_petsc  = K->pmatrix->matrix;
-		Vec Ml_petsc = Ml->pvector->vector;
-		Mat Mc_petsc = Mc->pmatrix->matrix;
-
-		/*Create D Matrix*/
-		#ifdef _HAVE_PETSC_
-		CreateDMatrix(&D_petsc,K_petsc);
-
-		/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
-		CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
-
-		/*Get previous solution u^n*/
-		GetSolutionFromInputsx(&ug,femmodel);
-		Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
-		delete ug;
-
-		/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
-		CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
-		delete uf;
-
-		/*Go solve lower order solution*/
-		SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
-		MatFree(&LHS);
-		VecFree(&RHS);
-
-		/*Richardson to calculate udot*/
-		RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
-
-		/*Serialize u and udot*/
-		IssmDouble* udot_serial = NULL;
-		IssmDouble* u_serial    = NULL;
-		IssmDouble* ml_serial   = NULL;
-		VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
-		VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
-		VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
-
-		/*Anti diffusive fluxes*/
-		Vec         Fbar            = NULL;
-		IssmDouble* Ri_minus_serial = NULL;
-		IssmDouble* Ri_plus_serial  = NULL;
-		CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,deltat);
-		CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
-
-		/*Clean up*/
-		MatFree(&D_petsc);
-		delete Mc;
-		xDelete<IssmDouble>(udot_serial);
-		xDelete<IssmDouble>(u_serial);
-		xDelete<IssmDouble>(ml_serial);
-		xDelete<IssmDouble>(Ri_plus_serial);
-		xDelete<IssmDouble>(Ri_minus_serial);
-
-		/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
-		Vec temp1 = NULL;
-		VecDuplicate(u,&temp1);
-		VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
-		VecAXPBY(udot,1.,1.,temp1);
-		VecAXPY(u,deltat,temp1);
-		VecFree(&Fbar);
-		VecFree(&udot);
-		VecFree(&temp1);
-
-		uf =new Vector<IssmDouble>(u);
-		VecFree(&u);
-
-		InputUpdateFromSolutionx(femmodel,uf); 
-		delete uf;
-
-		#else
-		_error_("PETSc needs to be installed");
-		#endif
-
-		delete Ml;
-		delete K;
-		delete analysis;
-
-	}
+void solutionsequence_fct(FemModel* femmodel){
+
+	/*intermediary: */
+	IssmDouble           theta,deltat,dmax;
+	int                  configuration_type;
+	Vector<IssmDouble>*  Ml = NULL;
+	Matrix<IssmDouble>*  K  = NULL;
+	Matrix<IssmDouble>*  Mc = NULL;
+	Vector<IssmDouble>*  ug = NULL;
+	Vector<IssmDouble>*  uf = NULL;
+
+	/*Create analysis*/
+	MasstransportAnalysis* analysis = new MasstransportAnalysis();
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->UpdateConstraintsx();
+	theta = 0.5;
+
+	/*Create lumped mass matrix*/
+	analysis->LumpedMassMatrix(&Ml,femmodel);
+	analysis->MassMatrix(&Mc,femmodel);
+	analysis->FctKMatrix(&K,NULL,femmodel);
+	delete analysis;
+
+	/*Convert matrices to PETSc matrices*/
+	Mat D_petsc  = NULL;
+	Mat LHS      = NULL;
+	Vec RHS      = NULL;
+	Vec u        = NULL;
+	Vec udot     = NULL;
+	Mat K_petsc  = K->pmatrix->matrix;
+	Vec Ml_petsc = Ml->pvector->vector;
+	Mat Mc_petsc = Mc->pmatrix->matrix;
+
+	/*Create D Matrix*/
+	#ifdef _HAVE_PETSC_
+	CreateDMatrix(&D_petsc,K_petsc);
+
+	/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
+	CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
+
+	/*Get previous solution u^n*/
+	GetSolutionFromInputsx(&ug,femmodel);
+	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
+	delete ug;
+
+	/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
+	CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
+	delete uf;
+
+	/*Go solve lower order solution*/
+	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
+	MatFree(&LHS);
+	VecFree(&RHS);
+
+	/*Richardson to calculate udot*/
+	RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
+	delete K;
+
+	/*Serialize u and udot*/
+	IssmDouble* udot_serial = NULL;
+	IssmDouble* u_serial    = NULL;
+	IssmDouble* ml_serial   = NULL;
+	VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
+	VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
+	VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
+
+	/*Anti diffusive fluxes*/
+	Vec         Fbar            = NULL;
+	IssmDouble *Ri_minus_serial = NULL;
+	IssmDouble *Ri_plus_serial  = NULL;
+	IssmDouble *ulmin           = NULL;
+	IssmDouble *ulmax           = NULL;
+	femmodel->GetInputLocalMinMaxOnNodesx(&ulmin,&ulmax,ThicknessEnum);
+	CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,ulmin,ulmax,deltat);
+	CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
+	xDelete<IssmDouble>(Ri_plus_serial);
+	xDelete<IssmDouble>(Ri_minus_serial);
+	xDelete<IssmDouble>(ulmin);
+	xDelete<IssmDouble>(ulmax);
+
+	/*Clean up*/
+	MatFree(&D_petsc);
+	delete Mc;
+	xDelete<IssmDouble>(udot_serial);
+	xDelete<IssmDouble>(u_serial);
+	xDelete<IssmDouble>(ml_serial);
+
+	/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
+	UpdateSolution(u,udot,Ml_petsc,Fbar,deltat);
+	uf =new Vector<IssmDouble>(u);
+	VecFree(&u);
+	VecFree(&Fbar);
+	VecFree(&udot);
+	delete Ml;
+
+	/*Update Element inputs*/
+	InputUpdateFromSolutionx(femmodel,uf); 
+	delete uf;
+
+	#else
+	_error_("PETSc needs to be installed");
+	#endif
+}