Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 20657)
+++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 20658)
@@ -385,5 +385,5 @@
 ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceSSAHO(IssmDouble kmax){/*{{{*/
 
-	const int numdof=NUMVERTICES*NDOF2;
+	int        numdof,numdof2,N;
 	IssmDouble penalty_offset;
 
@@ -394,14 +394,16 @@
 	parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
 
-	//Create elementary matrix: add penalty to 
-	Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
-	Ke->values[0*numdof+2]=-kmax*pow(10.,penalty_offset);
-	Ke->values[2*numdof+0]=-kmax*pow(10.,penalty_offset);
-	Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
-
-	Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
-	Ke->values[1*numdof+3]=-kmax*pow(10.,penalty_offset);
-	Ke->values[3*numdof+1]=-kmax*pow(10.,penalty_offset);
-	Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
+	/*Get number of dof for these two nodes*/
+	numdof =this->nodes[0]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+	numdof2=this->nodes[1]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+	N=NUMVERTICES*numdof;
+
+	/*Add penalty to Element matrix*/
+	for(int i=0;i<numdof;i++){
+		Ke->values[         i*N+i       ]=+kmax*pow(10.,penalty_offset);
+		Ke->values[         i*N+numdof+i]=-kmax*pow(10.,penalty_offset);
+		Ke->values[(numdof+i)*N+i       ]=-kmax*pow(10.,penalty_offset);
+		Ke->values[(numdof+i)*N+numdof+i]=+kmax*pow(10.,penalty_offset);
+	}
 
 	/*Clean up and return*/