Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16761)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16762)
@@ -947,8 +947,13 @@
 			return;
 		case L1L2ApproximationEnum:
-			InputUpdateFromSolutionHO(solution,element);
+			InputUpdateFromSolutionSSA(solution,element);
 			return;
-		case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
-			/*the elements around will create the solution*/
+		case SSAHOApproximationEnum:
+			InputUpdateFromSolutionSSAHO(solution,element);
+			return;
+		case HOFSApproximationEnum:
+			_error_("not implemented yet");
+		case SSAFSApproximationEnum:
+			_error_("not implemented yet");
 			return;
 		default:
@@ -975,10 +980,10 @@
 	switch(meshtype){
 		case Mesh2DhorizontalEnum:
-			element->GetInputListOnNodes(thickness,ThicknessEnum);
+			element->GetInputListOnVertices(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);
+			element->GetInputListOnVertices(surface,SurfaceEnum);
 			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
 			break;
@@ -996,5 +1001,5 @@
 			break;
 		case Mesh3DEnum:
-			if(!element->IsOnBed()) return;
+			if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
 			basalelement=element->SpawnBasalElement();
 			break;
@@ -1225,2 +1230,95 @@
 	xDelete<int>(cs_list);
 }/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i,meshtype;
+	IssmDouble  rho_ice,g;
+	int*        SSAdoflist = NULL;
+	int*        HOdoflist  = NULL;
+	IssmDouble* xyz_list   = NULL;
+
+	/*we have to add results of this element for HO and results from the element
+	 * at base for SSA, so we need to have the pointer toward the basal element*/
+	Element* basalelement=element->GetBasalElement();
+	if(basalelement->ObjectEnum()!=PentaEnum){
+		_error_("Coupling not supported for "<<EnumToStringx(basalelement->ObjectEnum()));
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+	int numdof2d = numnodes;/*FIXME: only works for Penta*/
+
+	/*Fetch dof list and allocate solution vectors*/
+	basalelement->GetDofList(&SSAdoflist,SSAApproximationEnum,GsetEnum);
+	element->GetDofList(     &HOdoflist, HOApproximationEnum, GsetEnum);
+	IssmDouble* HOvalues  = xNew<IssmDouble>(numdof);
+	IssmDouble* SSAvalues = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	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: */
+	for(i=0;i<numdof2d;i++){
+		HOvalues[i]  = solution[HOdoflist[i]];
+		SSAvalues[i] = solution[SSAdoflist[i]];
+	}
+	for(i=numdof2d;i<numdof;i++){
+		HOvalues[i]  = solution[HOdoflist[i]];
+		SSAvalues[i] = SSAvalues[i-numdof2d];
+	}
+
+	/*Transform solution in Cartesian Space*/
+	basalelement->TransformSolutionCoord(SSAvalues,XYEnum);
+	element->TransformSolutionCoord(HOvalues,XYEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numnodes;i++){
+		vx[i]=SSAvalues[i*2+0]+HOvalues[i*2+0];
+		vy[i]=SSAvalues[i*2+1]+HOvalues[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*/
+	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 
+	 * status, otherwise, we'll wipe them off: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+	element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddInput(VxEnum,vx,P1Enum);
+	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);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(SSAvalues);
+	xDelete<IssmDouble>(HOvalues);
+	xDelete<int>(SSAdoflist);
+	xDelete<int>(HOdoflist);
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16761)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16762)
@@ -25,4 +25,5 @@
 		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionSSAHO(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 16761)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16762)
@@ -53,4 +53,5 @@
 		virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
 		virtual void   TransformSolutionCoord(IssmDouble* values,int* transformenum_list)=0;
+		virtual Element* GetBasalElement(void)=0;
 		virtual void	GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
 		virtual void	GetDofListVelocity(int** pdoflist,int setenum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16761)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16762)
@@ -1029,5 +1029,5 @@
 /*}}}*/
 /*FUNCTION Penta::GetBasalElement{{{*/
-Penta* Penta::GetBasalElement(void){
+Element* Penta::GetBasalElement(void){
 
 	/*Output*/
@@ -5215,5 +5215,5 @@
 	
 	/*get basal element, needed for basal watercolumn*/
-	pentabase=this->GetBasalElement();
+	pentabase=(Penta*)this->GetBasalElement();
 	
 	this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
@@ -6988,5 +6988,5 @@
 
 	/*Find penta on bed as HO must be coupled to the dofs on the bed: */
-	Penta* pentabase=GetBasalElement();
+	Penta* pentabase=(Penta*)GetBasalElement();
 	Tria*  tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7183,5 +7183,5 @@
 
 	/*Find penta on bed as FS must be coupled to the dofs on the bed: */
-	Penta* pentabase=GetBasalElement();
+	Penta* pentabase=(Penta*)GetBasalElement();
 	Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7618,5 +7618,5 @@
 
 	/*Find penta on bed as this is a SSA elements: */
-	pentabase=GetBasalElement();
+	pentabase=(Penta*)GetBasalElement();
 	tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -7770,5 +7770,5 @@
 
 	/*Find penta on bed as this is a SSA elements: */
-	pentabase=GetBasalElement();
+	pentabase=(Penta*)GetBasalElement();
 	tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
@@ -9982,5 +9982,5 @@
 	/*OK, we have to add results of this element for HO 
 	 * and results from the penta at base for SSA. Now recover results*/
-	penta=GetBasalElement();
+	penta=(Penta*)GetBasalElement();
 
 	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
@@ -10075,5 +10075,5 @@
 	/*OK, we have to add results of this element for SSA 
 	 * and results from the penta at base for SSA. Now recover results*/
-	penta=GetBasalElement();
+	penta=(Penta*)GetBasalElement();
 
 	/*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16761)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16762)
@@ -85,4 +85,5 @@
 		void   Delta18oParameterization(void);
 		void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
+		Element* GetBasalElement(void);
 		void	 GetDofList(int** pdoflist,int approximation_enum,int setenum);
 		void	 GetDofListVelocity(int** pdoflist,int setenum);
@@ -228,5 +229,4 @@
 		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
 		Penta*         GetLowerElement(void);
-		Penta*         GetBasalElement(void);
 		void           InputChangeName(int input_enum, int enum_type_old);
 		void	         InputExtrude(int enum_type,int object_type);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16761)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16762)
@@ -87,4 +87,5 @@
 		void        FindParam(IssmDouble* pvalue,int paramenum);
 		int         FiniteElement(void);
+		Element*    GetBasalElement(void){_error_("not implemented yet");};
 		void        GetDofList(int** pdoflist,int approximation_enum,int setenum){_error_("not implemented yet");};
 		void        GetDofListVelocity(int** pdoflist,int setenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16761)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16762)
@@ -83,4 +83,5 @@
 		void        FindParam(IssmDouble* pvalue,int paramenum);
 		int         FiniteElement(void);
+		Element*    GetBasalElement(void){_error_("not implemented yet");};
 		void	      GetDofList(int** pdoflist,int approximation_enum,int setenum);
 		void	      GetDofListVelocity(int** pdoflist,int setenum);
