Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18234)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18235)
@@ -2881,4 +2881,5 @@
 	int         i,dim,epssize;
 	IssmDouble  r,rl,Jdet,viscosity,DU,DUl;
+	IssmDouble	normal[3];
 	IssmDouble *xyz_list = NULL;
 	IssmDouble *xyz_list_base = NULL;
@@ -2959,4 +2960,10 @@
 		element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
 		element->GetVerticesCoordinatesBase(&xyz_list_base);
+		element->NormalBase(&normal[0],xyz_list_base);
+
+		int         lsize;
+		IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
+		IssmDouble* C       = xNewZeroInit<IssmDouble>(dim*lnumdof);
+		IssmDouble* Cprime  = xNewZeroInit<IssmDouble>(dim*numdof);
 
 		delete gauss;
@@ -2966,12 +2973,17 @@
 
 			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-			this->GetCFS(C,element,dim,xyz_list_base,gauss);
-			this->GetCFSprime(Cprime,element,dim,xyz_list_base,gauss);
-			DUl = gauss->weight*Jdet*sqrt(rl);
-			TripleMultiply(C,1,lnumdof,1,
-						&DU,1,1,0,
-						Cprime,1,numdof,0,
+			this->GetCFS(C,element,dim,xyz_list,gauss);
+			this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
+			for(i=0;i<dim;i++)   Dlambda[i*epssize+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
+			TripleMultiply(C,dim,lnumdof,1,
+						Dlambda,dim,dim,0,
+						Cprime,dim,numdof,0,
 						CtCUzawa,1);
 		}
+		/*Delete base part*/
+		xDelete<IssmDouble>(xyz_list_base);
+		xDelete<IssmDouble>(Dlambda);
+		xDelete<IssmDouble>(C);
+		xDelete<IssmDouble>(Cprime);
 	}
 
@@ -3733,6 +3745,8 @@
 	int         i,dim;
 	IssmDouble  Jdet,r,pressure;
-	IssmDouble *xyz_list = NULL;
-	Gauss*      gauss    = NULL;
+	IssmDouble  bed_normal[3];
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
 
 	/*Get problem dimension*/
@@ -3757,4 +3771,5 @@
 	/*Get d and tau*/
 	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+	Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
 
 	gauss=element->NewGauss(5);
@@ -3773,4 +3788,33 @@
 			}
 		}
+	}
+
+	if(element->IsOnBase()){ 
+
+		IssmDouble   sigmann;
+		IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
+
+		element->GetVerticesCoordinatesBase(&xyz_list_base);
+		element->NormalBase(&bed_normal[0],xyz_list_base);
+
+		delete gauss;
+		Gauss* gauss=element->NewGaussBase(5);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+			gauss->GaussPoint(ig);
+
+			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+			element->NodalFunctionsVelocity(vbasis,gauss);
+			sigmann_input->GetInputValue(&sigmann, gauss);
+
+			for(i=0;i<numnodes;i++){
+				pe->values[i*dim+0] += - sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
+				pe->values[i*dim+1] += - sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
+				if(dim==3){
+					pe->values[i*dim+2]+= - sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
+				}
+			}
+		}
+		xDelete<IssmDouble>(xyz_list_base);
+		xDelete<IssmDouble>(vbasis);
 	}
 
@@ -4410,5 +4454,6 @@
 }/*}}}*/
 void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute C  matrix. C=[Cp1 Cp2 ...] where Cpi=phi_pi. 
+	/*Compute C  matrix. C=[Cp1 Cp2 ...] where:
+	 *     Cpi=[phi phi]. 
 	 */
 
@@ -4421,9 +4466,10 @@
 	/*Get nodal functions derivatives*/
 	IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
-	element->NodalFunctionsP1(basis,gauss);
+	element->NodalFunctions(basis,gauss);
 
 	/*Build B: */
 	for(int i=0;i<lnumnodes;i++){
-		C[i] = basis[i];
+		C[i*lnumnodes+0] = basis[i];
+		C[i*lnumnodes+1] = basis[i];
 	}
 
@@ -4432,36 +4478,50 @@
 }/*}}}*/
 void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute C'  matrix. C'=[C1' C2' C3'] 
-	 *			Ci' = [  dphi/dx   dphi/dy ]
+	/*	Compute C'  matrix. C'=[C1' C2' ...] 
+	 *			Ci' = [  phi  0  ]
+	 *			      [   0  phi ]
 	 *
 	 *	In 3d
-	 *     	   Ci=[ dh/dx   dh/dy    dh/dz  ]
+	 *			Ci' = [  phi  0   0  ]
+	 *			      [   0  phi  0  ]
+	 *			      [   0   0  phi ]
 	 *	where phi is the finiteelement function for node i.
 	 */
 
 	/*Fetch number of nodes for this finite element*/
-	int lnumnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(dim*lnumnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build C_prime: */
-	if(dim==2){
-		for(int i=0;i<lnumnodes;i++){
-			Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
-			Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
+	int vnumnodes = element->GetNumberOfNodes();
+	int vnumdof   = vnumnodes*dim;
+
+	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+	element->NodalFunctionsVelocity(vbasis,gauss);
+
+	/*Build B: */
+	if(dim==3){
+		for(int i=0;i<vnumnodes;i++){
+			Cprime[vnumdof*0+3*i+0] = vbasis[i];
+			Cprime[vnumdof*0+3*i+1] = 0.;
+			Cprime[vnumdof*0+3*i+2] = 0.;
+
+			Cprime[vnumdof*1+3*i+0] = 0.;
+			Cprime[vnumdof*1+3*i+1] = vbasis[i];
+			Cprime[vnumdof*1+3*i+2] = 0.;
+
+			Cprime[vnumdof*2+3*i+0] = 0.;
+			Cprime[vnumdof*2+3*i+1] = 0.;
+			Cprime[vnumdof*2+3*i+2] = vbasis[i];
 		}
 	}
 	else{
-		for(int i=0;i<lnumnodes;i++){
-			Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
-			Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
-			Cprime[dim*i+2] = dbasis[2*lnumnodes+i];
-		}
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(dbasis);
+		for(int i=0;i<vnumnodes;i++){
+			Cprime[vnumdof*0+2*i+0] = vbasis[i];
+			Cprime[vnumdof*0+2*i+1] = 0.;
+
+			Cprime[vnumdof*1+2*i+0] = 0.;
+			Cprime[vnumdof*1+2*i+1] = vbasis[i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(vbasis);
 }/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18234)
+++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18235)
@@ -12,4 +12,5 @@
 
 	parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRhopEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRholambdaEnum));
 }/*}}}*/
 void UzawaPressureAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -167,23 +168,59 @@
 void UzawaPressureAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
+	int        dim;
 	int        *doflist   = NULL;
+	IssmDouble rholambda,un;
+	IssmDouble bed_normal[3];
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes       = element->GetNumberOfNodes();
+	int numnodessigma  = element->NumberofNodes(P2Enum);
+	element->FindParam(&dim,DomainDimensionEnum);
 
 	/*Fetch dof list and allocate solution vector*/
 	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numnodes);
-	IssmDouble* pressure = 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);
 	element->GetInputListOnNodes(&pressure[0],PressureEnum);
 
+	/*Update pressure enum first*/
 	for(int i=0;i<numnodes;i++){
 		values[i]   = pressure[i] + solution[doflist[i]];
 	}
-
 	element->AddInput(PressureEnum,values,element->GetElementType());
 
+	/*Now compute sigmann if on base*/
+	element->GetInputListOnNodes(&sigmann[0],SigmaNNEnum);
+
+	if(element->IsOnBase()){ 
+
+		element->NormalBase(&bed_normal[0],xyz_list_tria);
+		element->FindParam(&rholambda,AugmentedLagrangianRholambdaEnum);
+
+		Gauss* gauss = element->NewGauss();
+		for(int i=0;i<numnodessigma;i++){
+			gauss->GaussNode(P2Enum,i);
+
+			vx_input->GetInputValue(&vx, gauss);
+			vy_input->GetInputValue(&vy, gauss);
+			un=bed_normal[0]*vx[i] + bed_normal[1]*vy[i];
+			if(dim==3){
+			   vz_input->GetInputValue(&vz, gauss);
+				un = un + bed_normal[2]*vz[i];
+			}
+			values[i] = sigmann[i] + rholambda*un;
+		}
+	}
+	element->AddInput(SigmaNNEnum,values,P2Enum));
+
 	/*Free ressources:*/
+	delete gauss;
 	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(sigmann);
 	xDelete<IssmDouble>(pressure);
 	xDelete<int>(doflist);