Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16773)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16774)
@@ -950,5 +950,6 @@
 			return;
 		case HOFSApproximationEnum:
-			_error_("not implemented yet");
+			InputUpdateFromSolutionHOFS(solution,element);
+			return;
 		case SSAFSApproximationEnum:
 			_error_("not implemented yet");
@@ -1006,5 +1007,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	element->TransformSolutionCoord(&values[0],&cs_list[0]);
+	element->TransformSolutionCoord(values,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
@@ -1127,4 +1128,101 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i;
+	IssmDouble  rho_ice,g,FSreconditioning;
+	int*        doflistHO  = NULL;
+	int*        doflistFSv = NULL;
+	int*        doflistFSp = NULL;
+
+	/*Only works with Penta for now*/
+	if(element->ObjectEnum()!=PentaEnum) _error_("Coupling not supported for "<<EnumToStringx(element->ObjectEnum()));
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes  = 6;
+	int numdofHO  = 6*2;
+	int numdofFSv = 6*3;
+	int numdofFSp = 6;
+
+	/*Fetch dof list and allocate solution vectors*/
+	element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+	element->GetDofList(&doflistHO, HOApproximationEnum, GsetEnum);
+	element->GetDofListPressure(&doflistFSp,GsetEnum);
+	IssmDouble* HOvalues  = xNew<IssmDouble>(numdofHO);
+	IssmDouble* FSvalues  = xNew<IssmDouble>(numdofFSv+numdofFSp);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vzHO      = xNew<IssmDouble>(numnodes);
+	IssmDouble* vzFS      = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(2*numnodes);
+	for(i=0;i<numnodes;i++) cs_list[i]          = XYZEnum;
+	for(i=0;i<numnodes;i++) cs_list[numnodes+i] = PressureEnum;
+
+	/*Use the dof list to index into the solution vector: */
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	for(i=0;i<numdofHO ;i++) HOvalues[i]=solution[doflistHO[i]];
+	for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
+	for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
+
+	/*Transform solution in Cartesian Space*/
+	element->TransformSolutionCoord(FSvalues,cs_list);
+	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]       = FSvalues[i*3+0]+HOvalues[i*2+0];
+		vy[i]       = FSvalues[i*3+1]+HOvalues[i*2+1];
+		vzFS[i]     = FSvalues[i*3+2];
+		pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get Vz and compute vel*/
+	element->GetInputListOnVertices(vzHO,VzHOEnum);
+	for(i=0;i<numnodes;i++){
+		vz[i] = vzHO[i]+vzFS[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);
+	element->InputChangeName(VzEnum,VzPicardEnum);
+	element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+	/*Add vx and vy as inputs to element: */
+	element->AddInput(VxEnum,vx,P1Enum);
+	element->AddInput(VyEnum,vy,P1Enum);
+	element->AddInput(VzEnum,vz,P1Enum);
+	element->AddInput(VzFSEnum,vzFS,P1Enum);
+	element->AddInput(VelEnum,vel,P1Enum);
+	element->AddInput(PressureEnum,pressure,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vzHO);
+	xDelete<IssmDouble>(vzFS);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(FSvalues);
+	xDelete<IssmDouble>(HOvalues);
+	xDelete<int>(doflistFSp);
+	xDelete<int>(doflistFSv);
+	xDelete<int>(doflistHO);
+	xDelete<int>(cs_list);
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16773)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16774)
@@ -24,4 +24,5 @@
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);