Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5922)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5923)
@@ -2078,46 +2078,32 @@
 ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){
 
-	int        i,j;
-	/* node data: */
-	const int  numdofp=NDOF2*NUMVERTICES;
-	const int  numdofs=NDOF4*NUMVERTICES;
-	int*       doflistp=NULL;
-	int*       doflists=NULL;
-
-	/*Matrices*/
-	double Kestokes_gg[24][24];
-	double Kepattyn_gg[12][12];
-	double Keps_gg[12][24]={0.0};
-	double Kesp_gg[24][12]={0.0};
-
-	/*Get Pattyn and Stokes doflist*/
-	GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); //Pattyn dof list
-	GetDofList(&doflists,StokesApproximationEnum,GsetEnum); //Stokes dof list
-
-	ISSMERROR("not supported yet");
-//	/*Get CreateKMatrixDiagnosticStokes*/
-//	Kestokes_gg=CreateKMatrixDiagnosticStokes;
-//	
-//	/*Modify it and plug in into Kgg*/
-//	for( i=0; i<numdofs; i++){
-//		for(j=0;j<NUMVERTICES; j++){
-//			Kesp_gg[i][(j-1)*NDOF2+0]+=Kestokes_gg[i][(j-1)*NDOF4+0];
-//			Kesp_gg[i][(j-1)*NDOF2+1]+=Kestokes_gg[i][(j-1)*NDOF4+1];
-//		}
-//	}
-//	MatSetValues(Kgg,numdofs,doflists,numdofp,doflistp,(const double*)Kesp_gg,ADD_VALUES);
-//
-//	/*Get CreateKMatrixDiagnosticPattyn*/
-//	Kepattyn_gg=CreateKMatrixDiagnosticPattyn;
-//
-//	/*Modify it and plug in into Kgg*/
-//	for( i=0; i<numdofp; i++){
-//		for(j=0;j<NUMVERTICES; j++){
-//			Kesp_gg[i][(j-1)*NDOF4+0]+=Kestokes_gg[i][(j-1)*NDOF2+0];
-//			Kesp_gg[i][(j-1)*NDOF4+1]+=Kestokes_gg[i][(j-1)*NDOF2+1];
-//		}
-//	}
-//	MatSetValues(Kgg,numdofp,doflistp,numdofs,doflists,(const double*)Keps_gg,ADD_VALUES);
-
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=NewElementMatrix(this->nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke2=NewElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1;
+	delete Ke2;
+	Ke1=CreateKMatrixDiagnosticPattyn();
+	Ke2=CreateKMatrixDiagnosticStokes();
+
+	/*Constants*/
+	const int    numdofp=2*NUMVERTICES;
+	const int    numdofs=4*NUMVERTICES;
+	const int    numdoftotal=(4+2)*NUMVERTICES;
+	int          i,j;
+
+	for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
+		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0];
+		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1];
+	}
+	for(i=0;i<numdofp;i++) for(j=0;j<NUMVERTICES;j++){
+		Ke->values[i*numdoftotal+numdofp+NDOF4*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
+		Ke->values[i*numdoftotal+numdofp+NDOF4*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
+	}
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
 }
 /*}}}*/
@@ -2529,5 +2515,5 @@
 	/*If on water or not Stokes, skip stiffness: */
 	inputs->GetParameterValue(&approximation,ApproximationEnum);
-	if(IsOnWater() || approximation!=StokesApproximationEnum) return NULL;
+	if(IsOnWater() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;
 	ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
 
@@ -2610,5 +2596,5 @@
 	/*If on water or not Stokes, skip stiffness: */
 	inputs->GetParameterValue(&approximation,ApproximationEnum);
-	if(IsOnWater() || IsOnShelf() || !IsOnBed() || approximation!=StokesApproximationEnum) return NULL;
+	if(IsOnWater() || IsOnShelf() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;
 	ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
 
@@ -4205,5 +4191,5 @@
 
 	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);