Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11340)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11341)
@@ -7437,4 +7437,6 @@
 		case StokesApproximationEnum:
 			return CreateJacobianDiagnosticStokes();
+		case NoneApproximationEnum:
+			return NULL;
 		default:
 			_error_("Approximation %s not supported yet",EnumToStringx(approximation));
@@ -7473,5 +7475,4 @@
 	/*Intermediaries */
 	int        i,j,ig;
-	int        approximation;
 	double     xyz_list[NUMVERTICES][3];
 	double     Jdet;
@@ -7489,5 +7490,4 @@
 
 	/*Retrieve all inputs and parameters*/
-	inputs->GetInputValue(&approximation,ApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
@@ -7535,30 +7535,29 @@
 ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
 
-	_error_("Not implemented yet");
 	/*Constants*/
-	const int    numdof=NDOF2*NUMVERTICES;
+	const int    numdof=NDOF4*NUMVERTICES;
 
 	/*Intermediaries */
 	int        i,j,ig;
-	int        approximation;
 	double     xyz_list[NUMVERTICES][3];
 	double     Jdet;
 	double     eps1dotdphii,eps1dotdphij;
 	double     eps2dotdphii,eps2dotdphij;
+	double     eps3dotdphii,eps3dotdphij;
 	double     mu_prime;
 	double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-	double     eps1[3],eps2[3];
+	double     eps1[3],eps2[3],eps3[3];
 	double     phi[NUMVERTICES];
 	double     dphi[3][NUMVERTICES];
 	GaussPenta *gauss=NULL;
 
-	/*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
-	ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+	/*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
+	ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
 
 	/*Retrieve all inputs and parameters*/
-	inputs->GetInputValue(&approximation,ApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -7573,7 +7572,7 @@
 		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
-		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
-		eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
-		eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
+		eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
+		eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
+		eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
 
 		for(i=0;i<6;i++){
@@ -7583,9 +7582,18 @@
 				eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
 				eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
-
-				Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
-				Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
-				Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
-				Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+				eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
+				eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
+
+				Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+				Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+				Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
+
+				Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+				Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+				Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
+
+				Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
+				Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
+				Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
 			}
 		}
@@ -7593,5 +7601,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
 
 	/*Clean up and return*/