Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18349)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18350)
@@ -111,6 +111,53 @@
 	*pdmax = dmax;
 	*pLHS  = LHS;
 }/*}}}*/
+void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
+	/*Create Left Hand side of Lower order solution
+	 *
+	 * RHS = [ML + (1 − theta) deltaT L^n] u^n
+	 *
+	 * where L = K + D
+	 *
+	 */
+
+	/*Intermediaries*/
+	Vec         Ku  = NULL;
+	Vec         Du  = NULL;
+	Vec         RHS = NULL;
+	int         dof;
+	IssmDouble  d;
+
+	/*Initialize vectors*/
+	VecDuplicate(u,&Ku);
+	VecDuplicate(u,&Du);
+	VecDuplicate(u,&RHS);
+
+	/*Create RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u*/
+	MatMult(K,u,Ku);
+	MatMult(D,u,Du);
+	VecPointwiseMult(RHS,Ml,u);
+	VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);
+	VecFree(&Ku);
+	VecFree(&Du);
+
+	/*Penalize Dirichlet boundary*/
+	for(int i=0;i<femmodel->constraints->Size();i++){
+		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		if(constraint->InAnalysis(configuration_type)){
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			d = d*dmax;
+			if(dof!=-1){
+				VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+			}
+		}
+	}
+	VecAssemblyBegin(RHS);
+	VecAssemblyEnd(  RHS);
+
+	/*Assign output pointer*/
+	*pRHS = RHS;
+
+}/*}}}*/
 #endif
 void solutionsequence_fct(FemModel* femmodel){
 
@@ -120,18 +167,15 @@
 	Matrix<IssmDouble>*  Mc = NULL;
 	Vector<IssmDouble>*  ug = NULL;
 	Vector<IssmDouble>*  uf = NULL;
-	Vector<IssmDouble>*  ys = NULL;
 
-	IssmDouble theta,deltat,dmax=0.;
-	int        dof,ncols,ncols2,ncols3,rstart,rend;
-	int        configuration_type,analysis_type;
-	double     d,mi;
+	IssmDouble theta,deltat,dmax;
+	int        dof,ncols,ncols2,rstart,rend;
+	int        configuration_type;
+	double     d;
 	int*       cols  = NULL;
 	int*       cols2 = NULL;
-	int*       cols3 = NULL;
 	double*    vals  = NULL;
 	double*    vals2 = NULL;
-	double*    vals3 = NULL;
 
 	/*Create analysis*/
 	MasstransportAnalysis* analysis = new MasstransportAnalysis();
@@ -139,7 +183,6 @@
 	/*Recover parameters: */
 	femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 	femmodel->UpdateConstraintsx();
 	theta = 0.5;
 
@@ -151,6 +194,8 @@
 	/*Convert matrices to PETSc matrices*/
 	Mat D_petsc  = NULL;
 	Mat LHS      = NULL;
+	Vec RHS      = NULL;
+	Vec u        = NULL;
 	Mat K_petsc  = K->pmatrix->matrix;
 	Vec Ml_petsc = Ml->pvector->vector;
 	Mat Mc_petsc = Mc->pmatrix->matrix;
@@ -166,35 +211,9 @@
 	GetSolutionFromInputsx(&ug,femmodel);
 	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
 	delete ug;
-	Vec u  = uf->pvector->vector;
-	Vec Ku = NULL;
-	Vec Du = NULL;
-	Vec RHS = NULL;
-	VecDuplicate(u,&Ku);
-	VecDuplicate(u,&Du);
-	VecDuplicate(u,&RHS);
-	MatMult(K_petsc,u,Ku);
-	MatMult(D_petsc,u,Du);
-	VecPointwiseMult(RHS,Ml_petsc,u);
-	VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);// RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u
-	VecFree(&Ku);
-	VecFree(&Du);
+	CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
 	delete uf;
 
-	/*Penalize Dirichlet boundary*/
-	for(int i=0;i<femmodel->constraints->Size();i++){
-		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-		if(constraint->InAnalysis(analysis_type)){
-			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-			d = d*dmax;
-			if(dof!=-1){
-				VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
-			}
-		}
-	}
-	VecAssemblyBegin(RHS);
-	VecAssemblyEnd(  RHS);
-
 	/*Go solve!*/
 	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
 	MatFree(&LHS);
@@ -324,9 +343,8 @@
 	uf =new Vector<IssmDouble>(u);
 	VecFree(&u);
 
-	Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
-	InputUpdateFromSolutionx(femmodel,ug); 
-	delete ug;  
+	InputUpdateFromSolutionx(femmodel,uf); 
+	delete uf;
 
 	#else
 	_error_("PETSc needs to be installed");