Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18280)
+++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18281)
@@ -31,4 +31,5 @@
 	if(iomodel->domaintype==Domain3DEnum) iomodel->FetchDataToInput(elements,VzEnum,0.);
 	iomodel->FetchDataToInput(elements,PressureEnum,0.);
+	iomodel->FetchDataToInput(elements,SigmaNNEnum,0.);
 }/*}}}*/
 void UzawaPressureAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -171,30 +172,27 @@
 	int        *doflist       = NULL;
 	IssmDouble rholambda,un,vx,vy,vz,sigmann;
-	IssmDouble bed_normal[3];
 	IssmDouble *xyz_list_base = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes       = element->GetNumberOfNodes();
-	//int numnodessigma  = element->NumberofNodes(P2Enum);
-	int numnodessigma;
+	int numnodes      = element->GetNumberOfNodes();
+	int numnodessigma = element->GetNumberOfNodes(P2Enum);
 
 	element->FindParam(&dim,DomainDimensionEnum);
-	if(dim==2) numnodessigma=3;
-	else numnodessigma=6;
 
 	/*Fetch dof list and allocate solution vector*/
 	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numnodes);
-	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
-	Input* vx_input=    element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input=    element->GetInput(VyEnum);       _assert_(vy_input);
-	Input* vz_input     = NULL;
-	if(dim==3){vz_input =element->GetInput(VzEnum);      _assert_(vz_input);}
-	Input* sigmann_input=element->GetInput(SigmaNNEnum); _assert_(vy_input);
+	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;
+	if(dim==3){vz_input       = element->GetInput(VzEnum);      _assert_(vz_input);}
 	element->GetInputListOnNodes(&pressure[0],PressureEnum);
 
 	/*Update pressure enum first*/
 	for(int i=0;i<numnodes;i++){
-		values[i]   = pressure[i] + solution[doflist[i]];
+		values[i]  = pressure[i] + solution[doflist[i]];
 	}
 	element->AddInput(PressureEnum,values,element->GetElementType());
@@ -202,33 +200,75 @@
 	/*Now compute sigmann if on base*/
 	if(element->IsOnBase()){ 
-
+		if(dim==3) _error_("not implemented yet");
+
+		int baselist[3];
+		int onbase=0;
+		IssmDouble  Jdet;
+		IssmDouble bed_normal[3];
+		IssmDouble  Jlambda[3][3]  = {0.0};
+		IssmDouble  Cuk[3]         = {0.0};
+		IssmDouble  deltalambda[3] = {0.0};
+		IssmDouble* vertexonbase  = xNew<IssmDouble>(numnodessigma);
+		Input* vertexonbase_input = element->GetInput(MeshVertexonbaseEnum); _assert_(vertexonbase_input);
+		Gauss* gauss = element->NewGauss();
+
+		IssmDouble* basis = xNewZeroInit<IssmDouble>(numnodessigma);
 		element->GetVerticesCoordinatesBase(&xyz_list_base);
 		element->NormalBase(&bed_normal[0],xyz_list_base);
 		element->FindParam(&rholambda,AugmentedLagrangianRholambdaEnum);
 
-		Gauss* gauss = element->NewGauss();
 		for(int i=0;i<numnodessigma;i++){
 			gauss->GaussNode(P2Enum,i);
-
-			sigmann_input->GetInputValue(&sigmann, gauss);
+			vertexonbase_input->GetInputValue(&vertexonbase[i], gauss);
+			if(vertexonbase[i]==1){ 
+				baselist[onbase]=i;
+				onbase += 1;
+			}
+		}
+		if(!onbase==3) _error_("basal nodes of element not found");
+
+		delete gauss;
+		gauss = element->NewGaussBase(3);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+			gauss->GaussPoint(ig);
+
+			/*Compute Jlambda*/
+			element->NodalFunctionsP2(basis,gauss);
+			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+			for(int i=0;i<3;i++){
+				for(int j=0;j<3;j++){
+					Jlambda[i][j] += Jdet*gauss->weight*basis[baselist[i]]*basis[baselist[j]];
+				}
+			}
+
+			/*Compute rho_lambd C u^k*/
 			vx_input->GetInputValue(&vx, gauss);
 			vy_input->GetInputValue(&vy, gauss);
 			un=bed_normal[0]*vx + bed_normal[1]*vy;
-			if(dim==3){
-			   vz_input->GetInputValue(&vz, gauss);
-				un = un + bed_normal[2]*vz;
-			}
-			values[i] = sigmann + rholambda*un;
-		}
+			for(int i=0;i<3;i++) Cuk[i] += - un*rholambda*Jdet*gauss->weight*basis[baselist[i]];
+		}
+
+		/*Now update sigmann*/
+		Matrix3x3Solve(&deltalambda[0],&Jlambda[0][0],&Cuk[0]);
 		delete gauss;
+		gauss = element->NewGauss();
+		for(int i=0;i<3;i++){
+			gauss->GaussNode(P2Enum,baselist[i]);
+			sigmann_input->GetInputValue(&sigmann, gauss);
+			valueslambda[baselist[i]] = sigmann + deltalambda[i];
+		}
+
+		delete gauss;
+		xDelete<IssmDouble>(vertexonbase);
 		xDelete<IssmDouble>(xyz_list_base);
-	}
-	element->AddInput(SigmaNNEnum,values,P2Enum);
+		xDelete<IssmDouble>(basis);
+	}
+	element->AddInput(SigmaNNEnum,valueslambda,P2Enum);
 
 	/*Free ressources:*/
 	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(valueslambda);
 	xDelete<IssmDouble>(pressure);
 	xDelete<int>(doflist);
-
 }/*}}}*/
 void UzawaPressureAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/