Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18206)
+++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18207)
@@ -92,11 +92,6 @@
 	int          dim;
 	IssmDouble   Jdet,r,divu;
-	IssmDouble   dvx[3],dvy[3],dvz[3];
 	IssmDouble   *xyz_list = NULL;
 	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element matrix and vectors*/
-	ElementVector* pe   = element->NewElementVector();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -105,7 +100,14 @@
 	element->GetVerticesCoordinates(&xyz_list);
 
+	/*Initialize Element matrix and vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dvx   = xNew<IssmDouble>(dim);
+	IssmDouble*    dvy   = xNew<IssmDouble>(dim);
+	IssmDouble*    dvz   = xNew<IssmDouble>(dim);
+
 	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
 	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input;
+	Input* vz_input = NULL;
 	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
 
@@ -122,8 +124,9 @@
 		}
 
+		divu=dvx[0]+dvy[1];
+		if (dim==3) divu=divu+dvz[2];
+
 		for(int i=0;i<numnodes;i++){
-			divu=dvx[0]+dvy[1];
-			if (dim==3) divu=divu+dvz[2];
-			pe->values[i]+=divu*r*Jdet*gauss->weight*basis[i];
+			pe->values[i] += - r * divu * Jdet * gauss->weight * basis[i];
 		}
 	}
@@ -133,4 +136,7 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dvx);
+	xDelete<IssmDouble>(dvy);
+	xDelete<IssmDouble>(dvz);
 	return pe;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 18206)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 18207)
@@ -13,5 +13,4 @@
 	Matrix<IssmDouble>*  Kff     = NULL;
 	Matrix<IssmDouble>*  Kfs     = NULL;
-	Vector<IssmDouble>*  ug_old  = NULL;
 	Vector<IssmDouble>*  ug      = NULL;
 	Vector<IssmDouble>*  uf      = NULL;
@@ -23,4 +22,6 @@
 	IssmDouble           eps_rel,r,theta; // 0<theta<.5   -> .15<theta<.45
 	int                  configuration_type,max_nonlinear_iterations;
+	bool                 vel_converged      = false;
+	bool                 pressure_converged = false;
 
 	/*Create analysis*/
@@ -40,8 +41,8 @@
 	/*Convergence criterion*/
 	int  count = 0;
-	GetSolutionFromInputsx(&ug,femmodel);
-	Vector<IssmDouble>* vx     = NULL;
-	Vector<IssmDouble>* vx_old = NULL;
-	GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+	Vector<IssmDouble>* vel     = NULL;
+	Vector<IssmDouble>* vel_old = NULL;
+	GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+	GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
 	while(true){
@@ -49,10 +50,11 @@
 
 		/*save pointer to old velocity*/
-		delete ug_old;ug_old=ug;
-		delete vx_old;vx_old=vx;
-		delete pug_old;pug_old=pug;
+		delete vel_old;vel_old=vel;
+		delete pug_old;pug_old=pug; 
+		pressure_converged=false; vel_converged=false;
 
 		/*Solve KU=F*/
 		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+		femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
 		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
@@ -63,8 +65,9 @@
 
 		/*Update solution*/
-		InputUpdateFromSolutionx(femmodel,ug); 
-		GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+		InputUpdateFromSolutionx(femmodel,ug); delete ug;
+		GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
 
 		femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum);
+		femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
 		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
@@ -72,25 +75,37 @@
 		Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
 		delete Kff; delete pf; delete df;
-		Mergesolutionfromftogx(&pug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+		Mergesolutionfromftogx(&pug, uf,ys,femmodel->nodes,femmodel->parameters); delete uf; delete ys;
 
 		/*Update solution*/
-		InputUpdateFromSolutionx(femmodel,pug); 
+		InputUpdateFromSolutionx(femmodel,pug); delete pug;
+		GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
 		/*Check for convergence*/
-		//Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
-		//IssmDouble ndu=dug->Norm(NORM_TWO);   delete dug;
-		//IssmDouble nu =ug_old->Norm(NORM_TWO);
-		Vector<IssmDouble>* dvx=vx_old->Duplicate(); vx_old->Copy(dvx); dvx->AYPX(vx,-1.0);
-		IssmDouble ndu=dvx->Norm(NORM_TWO);   delete dvx;
-		IssmDouble nu =vx_old->Norm(NORM_TWO);
+		Vector<IssmDouble>* dvel=vel_old->Duplicate(); vel_old->Copy(dvel); dvel->AYPX(vel,-1.0);
+		IssmDouble ndu=dvel->Norm(NORM_TWO);   delete dvel;
+		IssmDouble nu =vel_old->Norm(NORM_TWO);
 		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
 		if((ndu/nu)<eps_rel){
 			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
-			break;
+			vel_converged=true;
 		}
 		else{ 
 			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+			vel_converged=false;
+		}
+		Vector<IssmDouble>* dup=pug_old->Duplicate(); pug_old->Copy(dup); dup->AYPX(pug,-1.0);
+		IssmDouble ndp=dup->Norm(NORM_TWO);   delete dup;
+		IssmDouble np =pug_old->Norm(NORM_TWO);
+		if (xIsNan<IssmDouble>(ndp) || xIsNan<IssmDouble>(np)) _error_("convergence criterion is NaN!");
+		if((ndp/np)<eps_rel){
+			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(dp)/norm(p)" << ndp/np*100 << " < " << eps_rel*100 << " %\n");
+			pressure_converged=true;
+		}
+		else{ 
+			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(dp/)/norm(p)" << ndp/np*100 << " > " << eps_rel*100 << " %\n");
+			pressure_converged=false;
 		}
 
+		if(pressure_converged && vel_converged) break;
 		if(count>=max_nonlinear_iterations){
 			_printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
@@ -101,10 +116,8 @@
 	if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
 
-	delete ug;  
-	delete ug_old;  
 	delete pug;  
 	delete pug_old;  
-	delete vx;  
-	delete vx_old;  
+	delete vel;  
+	delete vel_old;  
 	delete stressanalysis;
 	delete pressureanalysis;