Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18348)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18349)
@@ -161,37 +161,5 @@
 
 	/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
-	MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
-	MatGetOwnershipRange(K_petsc,&rstart,&rend);
-	for(int row=rstart; row<rend; row++){
-		MatGetRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatGetRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		_assert_(ncols==ncols2);
-		for(int j=0; j<ncols; j++) {
-			_assert_(cols[j]==cols2[j]);
-			d = -theta*deltat*(vals[j] + vals2[j]);
-			if(cols[j]==row){
-				VecGetValues(Ml_petsc,1,(const int*)&cols[j],&mi);
-				d += mi;
-			}
-			if(fabs(d)>dmax) dmax = fabs(d);
-			MatSetValues(LHS,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
-		}
-		MatRestoreRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-	}
-
-	/*Penalize Dirichlet boundary*/
-	dmax = dmax * 1.e+3;
-	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);
-			if(dof!=-1){
-				MatSetValues(LHS,1,&dof,1,&dof,(const double*)&dmax,INSERT_VALUES);
-			}
-		}
-	}
-	MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+	CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
 
 	/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
@@ -265,4 +233,5 @@
 	VecDuplicate(u,&Ri_plus);
 	VecDuplicate(u,&Ri_minus);
+	MatGetOwnershipRange(K_petsc,&rstart,&rend);
 	for(int row=rstart; row<rend; row++){
 		double Pi_plus  = 0.;
