Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 5509)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 5510)
@@ -17,4 +17,5 @@
 	int segment_width;
 	int element;
+	int penpair_ids[2];
 	
 	/*DataSet*/
@@ -81,5 +82,5 @@
 	for (i=0;i<iomodel->numberofvertices;i++){
 		
-		if((iomodel->my_vertices[i]==1)){
+		if(iomodel->my_vertices[i]==1 && iomodel->singlenodetoelementconnectivity[i]!=0){
 
 			if ((iomodel->gridonbed[i]) && (iomodel->gridonicesheet[i]) && (iomodel->gridonstokes[i])){
@@ -97,4 +98,27 @@
 	xfree((void**)&iomodel->singlenodetoelementconnectivity);
 
+	/*Create Penpair for penalties: */
+	IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
+	
+	for(i=0;i<iomodel->numpenalties;i++){
+
+		if(iomodel->my_vertices[(int)iomodel->penalties[2*i+0]-1]){
+
+			/*In debugging mode, check that the second node is in the same cpu*/
+			ISSMASSERT(iomodel->my_vertices[(int)iomodel->penalties[2*i+1]-1]);
+
+			/*Get node ids*/
+			penpair_ids[0]=iomodel->nodecounter+(int)iomodel->penalties[2*i+0];
+			penpair_ids[1]=iomodel->nodecounter+(int)iomodel->penalties[2*i+1];
+
+			/*Create Load*/
+			loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],DiagnosticStokesAnalysisEnum));
+			count++;
+		}
+	}
+
+	/*free ressources: */
+	xfree((void**)&iomodel->penalties);
+
 	cleanup_and_return:
 
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5509)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5510)
@@ -2158,5 +2158,5 @@
 			}
 			if(!iomodel->pressure){
-				for(i=0;i<6;i++)nodeinputs[i]=i;
+				for(i=0;i<6;i++)nodeinputs[i]=0;
 				if(iomodel->qmu_analysis){
 					this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));
Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 5509)
+++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 5510)
@@ -44,4 +44,5 @@
 	ISSMASSERT(iomodel->singlenodetoelementconnectivity);
 	ISSMASSERT(index>=0 && index<iomodel->numberofvertices);
+	ISSMASSERT(id);
 
 	/*id: */
@@ -52,4 +53,5 @@
 	pengrid_node_id=iomodel->nodecounter+index+1;
 	pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
+	ISSMASSERT(pengrid_element_id);
 	pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
 
Index: /issm/trunk/src/c/objects/Loads/Penpair.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Penpair.cpp	(revision 5509)
+++ /issm/trunk/src/c/objects/Loads/Penpair.cpp	(revision 5510)
@@ -209,4 +209,7 @@
 		PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax);
 	}
+	else if (analysis_type==DiagnosticStokesAnalysisEnum){
+		PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax);
+	}
 	else{
 		ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
@@ -288,4 +291,57 @@
 }
 /*}}}1*/
+/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
+void  Penpair::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax){
+	
+	const int numgrids=2;
+	const int NDOF3=4;
+	const int numdof=numgrids*NDOF3;
+	int*      doflist=NULL;
+	int       numberofdofspernode;
+
+	double Ke[8][8]={0.0};
+	double penalty_offset;
+
+	/*pointers: */
+	Node** nodes=NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist);
+
+	/*recover pointers: */
+	nodes=(Node**)hnodes->deliverp();
+
+	/*recover parameters: */
+	parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
+
+	//Create elementary matrix: add penalty to 
+	Ke[0][0]=kmax*pow((double)10.0,penalty_offset);
+	Ke[0][4]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[4][0]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[4][4]=kmax*pow((double)10.0,penalty_offset);
+
+	Ke[1][1]=kmax*pow((double)10.0,penalty_offset);
+	Ke[1][5]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[5][1]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[5][5]=kmax*pow((double)10.0,penalty_offset);
+	
+	Ke[2][2]=kmax*pow((double)10.0,penalty_offset);
+	Ke[2][6]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[6][2]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[6][6]=kmax*pow((double)10.0,penalty_offset);
+
+	Ke[3][3]=kmax*pow((double)10.0,penalty_offset);
+	Ke[3][7]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[7][3]=-kmax*pow((double)10.0,penalty_offset);
+	Ke[7][7]=kmax*pow((double)10.0,penalty_offset);
+
+	/*Add Ke to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
+	
+	/*Free ressources:*/
+	xfree((void**)&doflist);
+
+}
+/*}}}1*/
 /*FUNCTION Penpair::GetDofList {{{1*/
 void  Penpair::GetDofList(int** pdoflist){
Index: /issm/trunk/src/c/objects/Loads/Penpair.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Penpair.h	(revision 5509)
+++ /issm/trunk/src/c/objects/Loads/Penpair.h	(revision 5510)
@@ -65,4 +65,5 @@
 			/*Penpair management: {{{1*/
 		void  PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax);
+		void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax);
 		void  GetDofList(int** pdoflist);
 		/*}}}*/