Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6999)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 7000)
@@ -719,16 +719,16 @@
 	const int    numdofm=NDOF2*NUMVERTICES2D;
 	const int    numdofs=NDOF4*NUMVERTICES;
-	const int    numdoftotal=numdofm+numdofs;
+	const int    numdoftotal=2*numdofm+numdofs;
 
 	/*Intermediaries */
 	int         i,j,ig;
 	double      Jdet;
-	double      viscosity; //viscosity
-	double      epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double      viscosity,stokesreconditioning; //viscosity
+	double      epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
 	double      xyz_list[NUMVERTICES][3];
-	double      B[4][numdofs];
+	double      B[4][numdofs+3];
 	double      Bprime[4][numdofm];
 	double      B2[3][numdofm];
-	double      Bprime2[3][numdofs];
+	double      Bprime2[3][numdofs+3];
 	double      D[4][4]={0.0};            // material matrix, simple scalar matrix.
 	double      D2[3][3]={0.0};            // material matrix, simple scalar matrix.
@@ -736,6 +736,6 @@
 	double      Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix 
 	double      Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix 
-	double      Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point.
-	double      Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point.
+	double      Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point.
+	double      Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point.
 	GaussPenta *gauss=NULL;
 	GaussTria  *gauss_tria=NULL;
@@ -754,4 +754,5 @@
 	/* Get node coordinates and dof list: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
@@ -777,8 +778,8 @@
 		D_scalar=2*viscosity*gauss->weight*Jdet;
 		for (i=0;i<3;i++) D[i][i]=D_scalar;
-		D[4][4]=(-1)*gauss->weight*Jdet;
+		D[3][3]=-gauss->weight*Jdet*stokesreconditioning;
 		for (i=0;i<3;i++) D2[i][i]=D_scalar;
 
-		TripleMultiply( &B[0][0],4,numdofs,1,
+		TripleMultiply( &B[0][0],4,numdofs+3,1,
 					&D[0][0],4,4,0,
 					&Bprime[0][0],4,numdofm,0,
@@ -787,5 +788,5 @@
 		TripleMultiply( &B2[0][0],3,numdofm,1,
 					&D2[0][0],3,3,0,
-					&Bprime2[0][0],3,numdofs,0,
+					&Bprime2[0][0],3,numdofs+3,0,
 					&Ke_gg_gaussian2[0][0],0);
 
@@ -794,6 +795,5 @@
 	} 
 	for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
-	for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[j][i];
-	_error_("coupling not finished yet");
+	for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
 
 	/*Clean-up and return*/
@@ -990,4 +990,5 @@
 	double      viscosity, oldviscosity, newviscosity, viscosity_overshoot;
 	double      epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double      epsilons[6]; //6 for stokes
 	double      xyz_list[NUMVERTICES][3];
 	double      B[3][numdof2d];
@@ -1040,6 +1041,6 @@
 		}
 		else if (approximation==MacAyealStokesApproximationEnum){
-			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-			matice->GetViscosity3dStokes(&newviscosity,&epsilon[0]);
+			this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+			matice->GetViscosity3dStokes(&newviscosity,&epsilons[0]);
 		}
 		else _error_("approximation %i (%s) not supported yet",approximation,EnumToString(approximation));
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 6999)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 7000)
@@ -134,5 +134,4 @@
 		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
 	}
-
 }
 /*}}}*/
@@ -245,7 +244,4 @@
 		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
 		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
 	}
 
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 6999)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 7000)
@@ -104,5 +104,5 @@
 			if (iomodel->vertices_type[io_index]==MacAyealStokesApproximationEnum && iomodel->gridonmacayeal[io_index]){
 				if(!iomodel->gridonbed[io_index]){
-					for(k=1;k<=gsize;k++) this->FreezeDof(k);
+					for(k=1;k<=2;k++) this->FreezeDof(k);
 				}
 			}