Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16699)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16700)
@@ -100,4 +100,115 @@
 }/*}}}*/
 void StressbalanceVerticalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
+
+	int          numnodes = element->GetNumberOfNodes();
+	int          numdof=numnodes*1;
+
+	int          i;
+	int          approximation;
+	int*         doflist  = NULL;
+	IssmDouble*  xyz_list = NULL;
+	IssmDouble   rho_ice,g;
+
+	IssmDouble*  values    = xNew<IssmDouble>(numdof);
+	IssmDouble*  vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vzSSA     = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vzHO      = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vzFS      = xNew<IssmDouble>(numnodes);
+	IssmDouble*  vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
+	IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
+
+	/*Get the approximation and do nothing if the element in FS or None*/
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+		return;
+	}
+
+	/*Get dof list and vertices coordinates: */
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Use the dof list to index into the solution vector vz: */
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+	for(i=0;i<numdof;i++){
+		vz[i]=values[i*1+0];
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get Vx and Vy*/
+	element->GetInputListOnNodes(&vx[0],VxEnum,0.0); //default is 0
+	element->GetInputListOnNodes(&vy[0],VyEnum,0.0); //default is 0
+
+	/*Do some modifications if we actually have a HOFS or SSAFS element*/
+	if(approximation==HOFSApproximationEnum){
+		Input* vzFS_input=element->GetInput(VzFSEnum);
+		if (vzFS_input){
+			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+			element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
+		}
+		else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
+		for(i=0;i<numnodes;i++){
+			vzHO[i]=vz[i];
+			vz[i]=vzHO[i]+vzFS[i];
+		}
+	}
+	else if(approximation==SSAFSApproximationEnum){
+		Input* vzFS_input=element->GetInput(VzFSEnum);
+		if (vzFS_input){
+			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+			element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
+		}
+		else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
+		for(i=0;i<numnodes;i++){
+			vzSSA[i]=vz[i];
+			vz[i]=vzSSA[i]+vzFS[i];
+		}
+	}
+
+	/*Now Compute vel*/
+	for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+
+	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
+	 *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
+	if(approximation!=HOFSApproximationEnum &&  approximation!=SSAFSApproximationEnum){
+		rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
+		g       = element->GetMaterialParameter(ConstantsGEnum);
+		element->GetInputListOnNodes(&surface[0],SurfaceEnum,0.);
+		for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+	}
+
+	/*Now, we have to move the previous Vz inputs to old 
+	 * status, otherwise, we'll wipe them off and add the new inputs: */
+	element->InputChangeName(VzEnum,VzPicardEnum);
+
+	if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
+		element->InputChangeName(PressureEnum,PressurePicardEnum);
+		element->AddInput(PressureEnum,pressure,P1Enum);
+	}
+	else if(approximation==HOFSApproximationEnum){
+		element->AddInput(VzHOEnum,vzHO,P1Enum);
+	}
+	else if(approximation==SSAFSApproximationEnum){
+		element->AddInput(VzSSAEnum,vzSSA,P1Enum);
+	}
+	element->AddInput(VzEnum,vz,P1Enum);
+	element->AddInput(VelEnum,vel,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(surface);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vzSSA);
+	xDelete<IssmDouble>(vzHO);
+	xDelete<IssmDouble>(vzFS);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+}/*}}}*/