Index: ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18653)
+++ ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18654)
@@ -80,24 +80,31 @@
 	ElementMatrix* Ke   = element->NewElementMatrix();
 	IssmDouble*    M    = xNew<IssmDouble>(numnodes);
 
+	IssmDouble connectivity;
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
+	int numvertices = element->GetNumberOfVertices();
 
-	Gauss* gauss = element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	//Gauss* gauss = element->NewGauss(5);
+	//for(int ig=gauss->begin();ig<gauss->end();ig++){
+	//	gauss->GaussPoint(ig);
 
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetM(M,element,gauss);
-		D_scalar=gauss->weight*Jdet;
-		TripleMultiply(M,1,numnodes,1,
-					&D_scalar,1,1,0,
-					M,1,numnodes,0,
-					&Ke->values[0],1);
+	//	element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+	//	this->GetM(M,element,gauss);
+	//	D_scalar=gauss->weight*Jdet;
+	//	//TripleMultiply(M,1,numnodes,1,
+	//	//			&D_scalar,1,1,0,
+	//	//			M,1,numnodes,0,
+	//	//			&Ke->values[0],1);
+
+	//}
+	for(int iv=0;iv<numvertices;iv++){
+		connectivity=(IssmDouble)element->VertexConnectivity(iv);
+		Ke->values[iv*numvertices+iv]=1./connectivity;
 	}
 
 	/*Clean up and return*/
-	delete gauss;
+	//delete gauss;
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(M);
 	return Ke;
@@ -201,7 +208,6 @@
 	IssmDouble* values        = xNew<IssmDouble>(numnodes);
 	IssmDouble* valueslambda  = xNewZeroInit<IssmDouble>(numnodessigma);
 	IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
-	Input* sigmann_input      = element->GetInput(SigmaNNEnum); _assert_(sigmann_input);
 	Input* vx_input           = element->GetInput(VxEnum);      _assert_(vx_input);
 	Input* vy_input           = element->GetInput(VyEnum);      _assert_(vy_input);
 	Input* vz_input           = NULL;
@@ -216,6 +222,7 @@
 
 	/*Now compute sigmann if on base*/
 	if(element->IsOnBase() && 0){ 
+		Input* sigmann_input      = element->GetInput(SigmaNNEnum); _assert_(sigmann_input);
 		if(dim==3) _error_("not implemented yet");
 
 		int baselist[3];
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 18653)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 18654)
@@ -38,9 +38,11 @@
 	femmodel->UpdateConstraintsx();
 
 	/*Convergence criterion*/
-	int  count = 0;
-	Vector<IssmDouble>* vel     = NULL;
-	Vector<IssmDouble>* vel_old = NULL;
+	int  count       = 0;
+	int  count_local = 0;
+	Vector<IssmDouble>* vel           = NULL;
+	Vector<IssmDouble>* vel_old       = NULL;
+	Vector<IssmDouble>* vel_old_local = NULL;
 	GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
 	GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
@@ -48,23 +50,44 @@
 		count++;
 
 		/*save pointer to old velocity*/
-		delete vel_old;vel_old=vel;
+		delete vel_old;vel_old=vel->Duplicate(); vel->Copy(vel_old);
 		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);
-		Reduceloadx(pf, Kfs, ys); delete Kfs;
-		Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
-		delete Kff; delete pf; delete df;
-		Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+		while(true){
+			count_local++;
+			delete vel_old_local;vel_old_local=vel;
+			/*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);
+			Reduceloadx(pf, Kfs, ys); delete Kfs;
+			Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
+			delete Kff; delete pf; delete df;
+			Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
 
-		/*Update solution*/
-		InputUpdateFromSolutionx(femmodel,ug); delete ug;
-		GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+			/*Update solution*/
+			InputUpdateFromSolutionx(femmodel,ug); delete ug;
+			GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+			/*Check for convergence*/
+			Vector<IssmDouble>* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0);
+			IssmDouble ndu=dvel_local->Norm(NORM_TWO);   delete dvel_local;
+			IssmDouble nu =vel_old_local->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;
+			}
+			else{ 
+				if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+			}
+			if(count_local>=max_nonlinear_iterations){
+				_printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+				break;
+			}
+		}
+		count_local=0;
 
 		femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum);
 		femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
@@ -80,9 +103,9 @@
 		GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
 
 		/*Check for convergence*/
-		Vector<IssmDouble>* dvel=vel_old->Duplicate(); vel_old->Copy(dvel); dvel->AYPX(vel,-1.0);
+		Vector<IssmDouble>* dvel=vel_old_local->Duplicate(); vel_old_local->Copy(dvel); dvel->AYPX(vel,-1.0);
 		IssmDouble ndu=dvel->Norm(NORM_TWO);   delete dvel;
-		IssmDouble nu =vel_old->Norm(NORM_TWO);
+		IssmDouble nu =vel_old_local->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");
@@ -118,6 +141,7 @@
 	delete pug_old;  
 	delete vel;  
 	delete vel_old;  
+	delete vel_old_local;  
 	delete stressanalysis;
 	delete pressureanalysis;
 }
