Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16779)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16780)
@@ -1060,8 +1060,20 @@
 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
-	int         i,meshtype;
-	IssmDouble  rho_ice,g;
+	int         i;
 	int*        doflist=NULL;
 	IssmDouble* xyz_list=NULL;
+
+	/*Deal with pressure first*/
+	int numvertices = element->GetNumberOfVertices();
+	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+	IssmDouble  rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  g         = element->GetMaterialParameter(ConstantsGEnum);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetInputListOnVertices(surface,SurfaceEnum);
+	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+	element->AddInput(PressureEnum,pressure,P1Enum);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(surface);
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -1076,6 +1088,4 @@
 	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
 	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
-	IssmDouble* surface   = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
@@ -1084,5 +1094,4 @@
 	/*Transform solution in Cartesian Space*/
 	element->TransformSolutionCoord(&values[0],XYEnum);
-	element->FindParam(&meshtype,MeshTypeEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -1099,12 +1108,4 @@
 	element->GetInputListOnNodes(&vz[0],VzEnum,0.);
 	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-
-	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
-	 *so the pressure is just the pressure at the bedrock: */
-	rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	g       = element->GetMaterialParameter(ConstantsGEnum);
-	element->GetVerticesCoordinates(&xyz_list);
-	element->GetInputListOnNodes(&surface[0],SurfaceEnum);
-	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 Vx and Vy inputs  to old 
@@ -1118,9 +1119,6 @@
 	element->AddInput(VyEnum,vy,P1Enum);
 	element->AddBasalInput(VelEnum,vel,P1Enum);
-	element->AddBasalInput(PressureEnum,pressure,P1Enum);
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(surface);
-	xDelete<IssmDouble>(pressure);
 	xDelete<IssmDouble>(vel);
 	xDelete<IssmDouble>(vz);