Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11337)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11338)
@@ -788,5 +788,5 @@
 #ifdef _HAVE_DIAGNOSTIC_
 		case DiagnosticHorizAnalysisEnum:
-			Ke=CreateJacobianDiagnosticPattyn();
+			Ke=CreateJacobianDiagnosticHoriz();
 			break;
 #endif
@@ -7424,7 +7424,93 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){
+
+	int approximation;
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+
+	switch(approximation){
+		case MacAyealApproximationEnum:
+			_error_("Not supported yet");
+		case PattynApproximationEnum:
+			return CreateJacobianDiagnosticPattyn();
+		case StokesApproximationEnum:
+			return CreateJacobianDiagnosticStokes();
+		default:
+			_error_("Approximation %s not supported yet",EnumToStringx(approximation));
+	}
+}
+/*}}}*/
 /*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
 ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
 
+	/*Constants*/
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	/*Intermediaries */
+	int        i,j,ig;
+	int        approximation;
+	double     xyz_list[NUMVERTICES][3];
+	double     Jdet;
+	double     eps1dotdphii,eps1dotdphij;
+	double     eps2dotdphii,eps2dotdphij;
+	double     mu_prime;
+	double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double     eps1[3],eps2[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();
+
+	/*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);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(5,5);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+
+		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];
+
+		for(i=0;i<6;i++){
+			for(j=0;j<6;j++){
+				eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
+				eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
+				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;
+			}
+		}
+	}
+
+	/*Transform Coordinate System*/
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
+
+	_error_("Not implemented yet");
 	/*Constants*/
 	const int    numdof=NDOF2*NUMVERTICES;
Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11337)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11338)
@@ -235,5 +235,7 @@
 		ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
 		ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
+		ElementMatrix* CreateJacobianDiagnosticHoriz(void);
 		ElementMatrix* CreateJacobianDiagnosticPattyn(void);
+		ElementMatrix* CreateJacobianDiagnosticStokes(void);
 		void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
 		void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11337)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11338)
@@ -78,5 +78,5 @@
 
 		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
-		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);//VecFree(&pJf);
+		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);VecFree(&pJf);
 		VecAXPY(uf,1.,duf);      VecFree(&duf);
 		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);