Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 16726)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 16727)
@@ -297,4 +297,4 @@
 	xDelete<IssmDouble>(phi);
 	xDelete<int>(doflist);
-	if(meshtype!=Mesh2DhorizontalEnum) delete element;
-}/*}}}*/
+	if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16726)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16727)
@@ -937,9 +937,15 @@
 			InputUpdateFromSolutionFS(solution,element);
 			return;
-		case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
-			InputUpdateFromSolutionHoriz(solution,element);
+		case SSAApproximationEnum: 
+			InputUpdateFromSolutionSSA(solution,element);
+			return;
+		case HOApproximationEnum: 
+			InputUpdateFromSolutionHO(solution,element);
+			return;
+		case SIAApproximationEnum:
+			InputUpdateFromSolutionHO(solution,element);
 			return;
 		case L1L2ApproximationEnum:
-			InputUpdateFromSolutionHoriz(solution,element);
+			InputUpdateFromSolutionHO(solution,element);
 			return;
 		case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
@@ -950,5 +956,104 @@
 	}
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i,meshtype;
+	IssmDouble  rho_ice,g;
+	int*        doflist=NULL;
+	IssmDouble* xyz_list=NULL;
+	Element*    basalelement=NULL;
+
+	/*Deal with pressure first*/
+	int numvertices = element->GetNumberOfVertices();
+	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+	g       =element->GetMaterialParameter(ConstantsGEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			element->GetInputListOnNodes(thickness,ThicknessEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+			break;
+		case Mesh3DEnum:
+			element->GetVerticesCoordinates(&xyz_list);
+			element->GetInputListOnNodes(surface,SurfaceEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+	element->AddInput(PressureEnum,pressure,P1Enum);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(thickness);
+	xDelete<IssmDouble>(surface);
+
+	/*Get basal element*/
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return;
+			basalelement=element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Fetch dof list and allocate solution vectors*/
+	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	basalelement->TransformSolutionCoord(&values[0],XYEnum);
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numnodes;i++){
+		vx[i]=values[i*2+0];
+		vy[i]=values[i*2+1];
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get Vz and compute vel*/
+	basalelement->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]);
+
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddBasalInput(VxEnum,vx,P1Enum);
+	element->AddBasalInput(VyEnum,vy,P1Enum);
+	element->AddBasalInput(VelEnum,vel,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+	if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,meshtype;
@@ -969,5 +1074,4 @@
 	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
 	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
-	IssmDouble* thickness = xNew<IssmDouble>(numnodes);
 	IssmDouble* surface   = xNew<IssmDouble>(numnodes);
 
@@ -997,16 +1101,7 @@
 	rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
 	g       = element->GetMaterialParameter(ConstantsGEnum);
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
-			for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
-			break;
-		case Mesh3DEnum:   
-			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]);
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
+	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 
@@ -1019,9 +1114,8 @@
 	element->AddInput(VxEnum,vx,P1Enum);
 	element->AddInput(VyEnum,vy,P1Enum);
-	element->AddInput(VelEnum,vel,P1Enum);
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddBasalInput(VelEnum,vel,P1Enum);
+	element->AddBasalInput(PressureEnum,pressure,P1Enum);
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(thickness);
 	xDelete<IssmDouble>(surface);
 	xDelete<IssmDouble>(pressure);
@@ -1033,5 +1127,4 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
-
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16726)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16727)
@@ -23,5 +23,6 @@
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16727)
@@ -60,4 +60,5 @@
 		virtual int    GetNumberOfNodesVelocity(void)=0;
 		virtual int    GetNumberOfNodesPressure(void)=0;
+		virtual int    GetNumberOfVertices(void)=0;
 		virtual void   GetNodesSidList(int* sidlist)=0;
 		virtual void   GetNodesLidList(int* sidlist)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16727)
@@ -139,9 +139,8 @@
 			Penta*     penta=NULL;
 
-			for(i=1;i<NUMVERTICES2D;i++){
+			for(i=0;i<NUMVERTICES2D;i++){
 				extrudedvalues[i]=values[i];
 				extrudedvalues[i+NUMVERTICES2D]=values[i];
 			}
-			this->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
 			penta=this;
 			for(;;){
@@ -1394,4 +1393,9 @@
 int Penta::GetNumberOfNodesVelocity(void){
 	return this->NumberofNodesVelocity();
+}
+/*}}}*/
+/*FUNCTION Penta::GetNumberOfVertices;{{{*/
+int Penta::GetNumberOfVertices(void){
+	return NUMVERTICES; 
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16727)
@@ -93,4 +93,5 @@
 		int    GetNumberOfNodesPressure(void);
 		int    GetNumberOfNodesVelocity(void);
+		int    GetNumberOfVertices(void);
 		IssmDouble GetMaterialParameter(int enum_in);
 		void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16727)
@@ -98,4 +98,5 @@
 		int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
 		int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
+		int         GetNumberOfVertices(void){_error_("not implemented yet");};
 		void        GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");};
 		int         Sid(){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16727)
@@ -1367,4 +1367,9 @@
 int Tria::GetNumberOfNodes(void){
 	return this->NumberofNodes();
+}
+/*}}}*/
+/*FUNCTION Tria::GetNumberOfVertices;{{{*/
+int Tria::GetNumberOfVertices(void){
+	return NUMVERTICES;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16726)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16727)
@@ -91,4 +91,5 @@
 		int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
 		int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
+		int         GetNumberOfVertices(void);
 		int         Sid();
 		bool        IsOnBed();