Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15452)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15453)
@@ -3119,5 +3119,5 @@
 	/*Initialize Element matrix, vectors and Gaussian points*/
 	ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); //Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)
-	IssmDouble*    dphi   = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -3134,5 +3134,5 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsDerivatives(dphi,&xyz_list[0][0],gauss);
+		GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
 
 		thickness_input->GetInputValue(&thickness, gauss);
@@ -3144,8 +3144,8 @@
 		for(i=0;i<numnodes;i++){
 			for(j=0;j<numnodes;j++){
-				eps1dotdphii=eps1[0]*dphi[0*numnodes+i]+eps1[1]*dphi[1*numnodes+i];
-				eps1dotdphij=eps1[0]*dphi[0*numnodes+j]+eps1[1]*dphi[1*numnodes+j];
-				eps2dotdphii=eps2[0]*dphi[0*numnodes+i]+eps2[1]*dphi[1*numnodes+i];
-				eps2dotdphij=eps2[0]*dphi[0*numnodes+j]+eps2[1]*dphi[1*numnodes+j];
+				eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
+				eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
+				eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
+				eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
 
 				Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
@@ -3158,9 +3158,9 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+	TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
-	xDelete<IssmDouble>(dphi);
 	return Ke;
 }
@@ -5214,18 +5214,18 @@
 
 	/*Constants*/
-	const int    numdof=NDOF2*NUMVERTICES;
+	const int numnodes = this->GetNumberOfNodes();
+	const int numdof   = NDOF2*numnodes;
 
 	/*Intermediaries */
-	int        i,j;
-	bool       incomplete_adjoint;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdet,thickness;
-	IssmDouble eps1dotdphii,eps1dotdphij;
-	IssmDouble eps2dotdphii,eps2dotdphij;
-	IssmDouble mu_prime;
-	IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/
-	IssmDouble eps1[2],eps2[2];
-	IssmDouble dphi[2][NUMVERTICES];
-	GaussTria *gauss=NULL;
+	int         i,j;
+	bool        incomplete_adjoint;
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  Jdet,thickness;
+	IssmDouble  eps1dotdphii,eps1dotdphij;
+	IssmDouble  eps2dotdphii,eps2dotdphij;
+	IssmDouble  mu_prime;
+	IssmDouble  epsilon[3];/* epsilon=[exx,eyy,exy];*/
+	IssmDouble  eps1[2],eps2[2];
+	GaussTria  *gauss=NULL;
 
 	/*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
@@ -5240,4 +5240,7 @@
 	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
 
+	/*Allocate dbasis*/
+	IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
+
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
@@ -5247,5 +5250,5 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+		GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
 
 		thickness_input->GetInputValue(&thickness, gauss);
@@ -5255,15 +5258,15 @@
 		eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
 
-		for(i=0;i<3;i++){
-			for(j=0;j<3;j++){
-				eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
-				eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
-				eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
-				eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
-
-				Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
-				Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
-				Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
-				Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
+		for(i=0;i<numnodes;i++){
+			for(j=0;j<numnodes;j++){
+				eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
+				eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
+				eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
+				eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
+
+				Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
+				Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
+				Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
+				Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
 			}
 		}
@@ -5271,8 +5274,9 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+	TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
 	delete gauss;
+	xDelete<IssmDouble>(dbasis);
 	//Ke->Transpose();
 	return Ke;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 15452)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 15453)
@@ -73,6 +73,6 @@
 
 		pJf=pf->Duplicate();
-		Kff->MatMult(uf,pJf);// delete Kff);
-		pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     //delete pf);
+		Kff->MatMult(uf,pJf);
+		pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);
 
 		femmodel->CreateJacobianMatrixx(&Jff,kmax);