Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 24713)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 24714)
@@ -8,4 +8,5 @@
 #include "../modules/modules.h"
 #include "../analyses/analyses.h"
+#define USEPENALTYMETHOD false
 
 #ifdef _HAVE_PETSC_
@@ -94,19 +95,45 @@
 		MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
 	}
-
-	/*Broadcast max(dmax)*/
-	IssmDouble dmax_all;
-	ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	dmax = dmax_all;
-
-	/*Penalize Dirichlet boundary*/
-	dmax = dmax * 1.e+3;
-	for(int i=0;i<femmodel->constraints->Size();i++){
-		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-		constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-		if(dof!=-1){
-			MatSetValue(LHS,dof,dof,dmax,INSERT_VALUES);
-		}
+	MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+
+	/*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
+	if(USEPENALTYMETHOD){
+		/*Option 1: Penalty method*/
+
+		/*Broadcast max(dmax)*/
+		IssmDouble dmax_all;
+		ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
+		ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+		dmax = dmax_all;
+
+		dmax = dmax * 1.e+3;
+		for(int i=0;i<femmodel->constraints->Size();i++){
+			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			if(dof!=-1){
+				MatSetValue(LHS,dof,dof,dmax,INSERT_VALUES);
+			}
+		}
+	}
+	else{
+		/*Option 2: zero stiffness matrix, 1 one diagonal*/
+		int numrows = femmodel->constraints->Size();
+		int* rows = xNew<int>(numrows);
+		IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
+		numrows = 0;
+
+		dmax = dmax * 1.e+3;
+		for(int i=0;i<femmodel->constraints->Size();i++){
+			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			if(dof!=-1){
+				rows[numrows] = dof;
+				numrows++;
+			}
+		}
+		MatZeroRows(LHS,numrows,rows,1.,NULL,NULL);
+		xDelete<int>(rows);
+		xDelete<IssmDouble>(rows_spc);
 	}
 	MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
@@ -146,12 +173,37 @@
 	VecFree(&Du);
 
-	/*Penalize Dirichlet boundary*/
-	for(int i=0;i<femmodel->constraints->Size();i++){
-		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-		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);
+
+	/*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
+	if(USEPENALTYMETHOD){
+		/*Option 1: Penalty method*/
+		for(int i=0;i<femmodel->constraints->Size();i++){
+			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			d = d*dmax;
+			if(dof!=-1){
+				VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+			}
+		}
+	}
+	else{
+		/*Option 2: zero stiffness matrix, 1 one diagonal*/
+		int  numrows = femmodel->constraints->Size();
+		int* rows = xNew<int>(numrows);
+		IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
+		numrows = 0;
+		for(int i=0;i<femmodel->constraints->Size();i++){
+			Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			if(dof!=-1){
+				rows[numrows] = dof;
+				rows_spc[numrows] = d;
+				numrows++;
+			}
+		}
+		VecSetValues(RHS,numrows,rows,rows_spc,INSERT_VALUES);
+		xDelete<int>(rows);
+		xDelete<IssmDouble>(rows_spc);
 	}
 	VecAssemblyBegin(RHS);