Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5595)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5596)
@@ -89,4 +89,5 @@
 	NoneApproximationEnum,
 	PattynApproximationEnum,
+	PattynStokesApproximationEnum,
 	StokesApproximationEnum,
 	/*}}}*/
Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5595)
+++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5596)
@@ -79,4 +79,5 @@
 		case NoneApproximationEnum : return "NoneApproximation";
 		case PattynApproximationEnum : return "PattynApproximation";
+		case PattynStokesApproximationEnum : return "PattynStokesApproximation";
 		case StokesApproximationEnum : return "StokesApproximation";
 		case ObjectEnum : return "Object";
Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5595)
+++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5596)
@@ -77,4 +77,5 @@
 	else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
 	else if (strcmp(name,"PattynApproximation")==0) return PattynApproximationEnum;
+	else if (strcmp(name,"PattynStokesApproximation")==0) return PattynStokesApproximationEnum;
 	else if (strcmp(name,"StokesApproximation")==0) return StokesApproximationEnum;
 	else if (strcmp(name,"Object")==0) return ObjectEnum;
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 5595)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 5596)
@@ -83,4 +83,44 @@
 					else ISSMERROR("if vertices_type is MacAyealPattyn, you shoud have gridonpattyn or gridonmacayeal");
 			}
+			/*Also add spcs of coupling: zero at the border pattyn/stokes for the appropriate dofs*/
+			else if ((int)iomodel->vertices_type[i]==PattynStokesApproximationEnum){
+				/*If grionmacayeal, spc pattyn dofs: 3 & 4*/
+					if ((int)iomodel->gridonpattyn[i]){
+						constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+						constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+						constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+						if ((int)iomodel->spcvelocity[6*i+0]){
+							constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+							count++;
+						}
+						if ((int)iomodel->spcvelocity[6*i+1]){
+							constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+							count++;
+						}
+
+					}
+					else if ((int)iomodel->gridonstokes[i]){
+						constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+						constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+						if ((int)iomodel->spcvelocity[6*i+0]){
+							constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,*(iomodel->spcvelocity+6*i+3)/iomodel->yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+							count++;
+						}
+						if ((int)iomodel->spcvelocity[6*i+1]){
+							constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,*(iomodel->spcvelocity+6*i+4)/iomodel->yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+							count++;
+						}
+						if ((int)iomodel->spcvelocity[6*i+2]){
+							constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,*(iomodel->spcvelocity+6*i+5)/iomodel->yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+							count++;
+						}
+					}
+					else ISSMERROR("if vertices_type is PattynStokes, you shoud have gridonpattyn or gridonstokes");
+			}
 			/*Now add the regular spcs*/
 			else{
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 5595)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 5596)
@@ -76,4 +76,10 @@
 			count++;
 			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
+			count++;
+		}
+		else if ((int)*(iomodel->elements_type+element)==(PattynStokesApproximationEnum)){
+			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
+			count++;
+			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
 			count++;
 		}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp	(revision 5595)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp	(revision 5596)
@@ -41,4 +41,14 @@
 			doftype[3]=PattynApproximationEnum;
 		}
+		else if (vertices_type[0]==PattynStokesApproximationEnum){
+			numdofs=6;
+			doftype=(int*)xmalloc(numdofs*sizeof(int));
+			doftype[0]=PattynApproximationEnum;
+			doftype[1]=PattynApproximationEnum;
+			doftype[2]=StokesApproximationEnum;
+			doftype[3]=StokesApproximationEnum;
+			doftype[4]=StokesApproximationEnum;
+			doftype[5]=StokesApproximationEnum;
+		}
 		else ISSMERROR("Approximationtype %i (%s) not implemented yet for DiagnosticHoriz",(int)*vertices_type,EnumToString((int)*vertices_type));
 	}
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5595)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5596)
@@ -795,4 +795,9 @@
 			CreateKMatrixCouplingMacAyealPattyn( Kgg);
 		}
+		else if(approximation==PattynStokesApproximationEnum){
+			CreateKMatrixDiagnosticPattyn( Kgg);
+			CreateKMatrixDiagnosticStokes( Kgg);
+			//CreateKMatrixCouplingPattynStokes( Kgg);
+		}
 		else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
 	}
@@ -876,4 +881,9 @@
 			CreatePVectorDiagnosticMacAyeal( pg);
 			CreatePVectorDiagnosticPattyn( pg);
+		}
+		else if(approximation==PattynStokesApproximationEnum){
+			CreatePVectorDiagnosticPattyn( pg);
+			CreatePVectorDiagnosticStokes( pg);
+			//CreatePVectorDiagnosticPattynStokes( pg);
 		}
 		else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
@@ -3006,5 +3016,5 @@
 	/* Get node coordinates and dof list: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-	GetDofList(&doflist);
+	GetDofList(&doflist,StokesApproximationEnum);
 
 	/* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -4182,5 +4192,5 @@
 	/* Get node coordinates and dof list: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-	GetDofList(&doflist);
+	GetDofList(&doflist,StokesApproximationEnum);
 
 	/* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5595)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5596)
@@ -705,5 +705,5 @@
 
 	/* Get dof list and node coordinates: */
-	GetDofList(&doflist);
+	GetDofList(&doflist,StokesApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	
Index: /issm/trunk/src/c/objects/Loads/Penpair.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Penpair.cpp	(revision 5595)
+++ /issm/trunk/src/c/objects/Loads/Penpair.cpp	(revision 5596)
@@ -231,4 +231,12 @@
 				switch(approximation1){
 					case StokesApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break;
+					case NoneApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break;
+					default: ISSMERROR("not supported yet");
+				}
+				break;
+			case NoneApproximationEnum:
+				switch(approximation1){
+					case StokesApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break;
+					case NoneApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break;
 					default: ISSMERROR("not supported yet");
 				}
@@ -351,13 +359,13 @@
 	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);
+	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: */