Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16692)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16693)
@@ -930,4 +930,93 @@
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	switch(approximation){
+		case FSApproximationEnum: case NoneApproximationEnum:
+			//InputUpdateFromSolutionFS(solution,element);
+			return;
+		case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
+			InputUpdateFromSolutionHoriz(solution,element);
+			return;
+		case L1L2ApproximationEnum:
+			//InputUpdateFromSolutionHoriz(solution,element);
+			return;
+		case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
+			/*the elements around will create the solution*/
+			return;
+		default:
+			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
+	}
 }/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
+
+	int        i;
+	IssmDouble rho_ice,g;
+	int*       doflist=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Fetch dof list and allocate solution vectors*/
+	element->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);
+	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+	IssmDouble* thickness = 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*/
+	element->TransformSolutionCoord(&values[0],XYEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numnodes;i++){
+		vx[i]=values[i*NDOF2+0];
+		vy[i]=values[i*NDOF2+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: */
+	Matpar* matpar=element->GetMatparPointer();
+	rho_ice=matpar->GetRhoIce();
+	g=matpar->GetG();
+	element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
+	for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[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(PressureEnum,PressurePicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddInput(VxEnum,vx,P1Enum);
+	element->AddInput(VyEnum,vy,P1Enum);
+	element->AddInput(VelEnum,vel,P1Enum);
+	element->AddInput(PressureEnum,pressure,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(thickness);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<int>(doflist);
+
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16692)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16693)
@@ -19,7 +19,8 @@
 		void CreateLoads(Loads* loads, IoModel* iomodel);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
 		void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
+		void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolution(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 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16693)
@@ -21,4 +21,5 @@
 class Vertices;
 class Materials;
+class Matpar;
 class Input;
 class Gauss;
@@ -38,4 +39,5 @@
 		virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
 		virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
+		virtual void   AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
 		virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
@@ -45,4 +47,6 @@
 		virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;
 		virtual int    FiniteElement(void)=0;
+		virtual Matpar* GetMatparPointer(void)=0;
+		virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
 		virtual void	GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
 		virtual void	GetDofListVelocity(int** pdoflist,int setenum)=0;
@@ -58,4 +62,6 @@
 		virtual bool   IsNodeOnShelfFromFlags(IssmDouble* flags)=0; 
 		virtual bool   IsOnBed()=0;
+		virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype)=0;
+		virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
 		virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype)=0;
 		virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
@@ -69,4 +75,5 @@
 		virtual IssmDouble SurfaceArea(void)=0;
 		virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;
+		virtual void   InputChangeName(int enum_type,int enum_type_old)=0;
 		virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
 		virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16693)
@@ -1815,4 +1815,11 @@
 }
 /*}}}*/
+/*FUNCTION Penta::InputChangeName{{{*/
+void  Penta::InputChangeName(int new_enum,int original_enum){
+
+	/*Call inputs method*/
+	this->inputs->ChangeEnum(original_enum,new_enum);
+}
+/*}}}*/
 /*FUNCTION Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
 void Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){
@@ -3227,4 +3234,11 @@
 	return dt;
 }/*}}}*/
+/*FUNCTION Penta::TransformSolutionCoord{{{*/
+void Penta::TransformSolutionCoord(IssmDouble* values,int transformenum){
+
+	::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
 /*FUNCTION Penta::Update {{{*/
 void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ 
@@ -6286,5 +6300,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
+	::TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
 
 	/*fill in all arrays: */
@@ -9799,5 +9813,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
+	::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
@@ -9898,6 +9912,6 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
+	::TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
+	::TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -9991,6 +10005,6 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
+	::TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
+	::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -10073,5 +10087,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
+	::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
@@ -10155,5 +10169,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
+	::TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
@@ -10245,6 +10259,6 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
+	::TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
+	::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -10506,5 +10520,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
+	::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16693)
@@ -91,4 +91,5 @@
 		void   GetNodesLidList(int* lidlist);
 		int    GetNumberOfNodes(void);
+		Matpar* GetMatparPointer(void){_error_("not implemented yet");};
 		void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
 		IssmDouble GetZcoord(GaussPenta* gauss);
@@ -179,4 +180,5 @@
 		/*}}}*/
 		/*Penta specific routines:{{{*/
+		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
 		void	         BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
 		ElementMatrix* CreateBasalMassMatrix(void);
@@ -218,4 +220,5 @@
 		Penta*         GetLowerElement(void);
 		Penta*         GetBasalElement(void);
+		void           InputChangeName(int input_enum, int enum_type_old);
 		void	         InputExtrude(int enum_type,int object_type);
 		void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
@@ -235,4 +238,5 @@
 		Tria*	         SpawnTria(int location);
 		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
+		void           TransformSolutionCoord(IssmDouble* values,int transformenum);
 
 		#ifdef _HAVE_STRESSBALANCE_
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16693)
@@ -66,4 +66,5 @@
 		/*}}}*/
 		/*Element virtual functions definitions: {{{*/
+		void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
 		void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
@@ -84,4 +85,7 @@
 		void        GetDofListVelocity(int** pdoflist,int setenum){_error_("not implemented yet");};
 		void        GetDofListPressure(int** pdoflist,int setenum){_error_("not implemented yet");};
+		void        GetInputListOnNodes(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
+		void        GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");};
+		Matpar*     GetMatparPointer(void){_error_("not implemented yet");};
 		int         GetNodeIndex(Node* node){_error_("not implemented yet");};
 		void        GetNodesSidList(int* sidlist){_error_("not implemented yet");};
@@ -89,4 +93,5 @@
 		int         GetNumberOfNodes(void){_error_("not implemented yet");};
 		int         Sid(){_error_("not implemented yet");};
+		void        InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");};
 		bool        IsOnBed(){_error_("not implemented yet");};
 		bool        IsFloating(){_error_("not implemented yet");};
@@ -117,5 +122,4 @@
 		void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
 		void    ComputeEPLThickness(void){_error_("not implemented yet");};
-		void    UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");};
 		#endif
 		void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
@@ -135,6 +139,8 @@
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
 		IssmDouble  TimeAdapt(){_error_("not implemented yet");};
+		void        TransformSolutionCoord(IssmDouble* values,int transformenum){_error_("not implemented yet");};
 		void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
 		void UpdateConstraintsExtrudeFromTop(){_error_("not implemented");};
+		void UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");};
 
 #ifdef _HAVE_RESPONSES_
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16693)
@@ -164,4 +164,11 @@
 	*pd_nz=d_nz;
 	*po_nz=o_nz;
+}
+/*}}}*/
+/*FUNCTION Tria::AddInput{{{*/
+void  Tria::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+	/*Call inputs method*/
+	this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
@@ -1490,4 +1497,9 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetMatparPointer(void) {{{*/
+Matpar* Tria::GetMatparPointer(void){
+	return this->matpar;
+}
+/*}}}*/
 /*FUNCTION Tria::GetVertexPidList {{{*/
 void  Tria::GetVertexPidList(int* doflist){
@@ -1598,4 +1610,11 @@
 		inputs->DuplicateInput(original_enum,new_enum);
 	}
+}
+/*}}}*/
+/*FUNCTION Tria::InputChangeName{{{*/
+void  Tria::InputChangeName(int new_enum,int original_enum){
+
+	/*Call inputs method*/
+	this->inputs->ChangeEnum(original_enum,new_enum);
 }
 /*}}}*/
@@ -2765,4 +2784,11 @@
 }
 /*}}}*/
+/*FUNCTION Tria::TransformSolutionCoord{{{*/
+void Tria::TransformSolutionCoord(IssmDouble* values,int transformenum){
+
+	::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
 /*FUNCTION Tria::Update{{{*/
 void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){
@@ -4339,5 +4365,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
+	::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -4421,5 +4447,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
+	::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
@@ -4488,5 +4514,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
+	::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -6258,5 +6284,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
+	::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16692)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16693)
@@ -191,4 +191,5 @@
 		/*}}}*/
 		/*Tria specific routines:{{{*/
+		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		ElementMatrix* CreateKMatrix(void);
 		ElementMatrix* CreateKMatrixBalancethickness(void);
@@ -237,4 +238,5 @@
 		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
 		void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
+      Matpar*        GetMatparPointer(void);
 		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
 		Input*         GetInput(int inputenum);
@@ -250,4 +252,5 @@
 		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
+		void	         InputChangeName(int enum_type,int enum_type_old);
 		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
 		void	         InputUpdateFromSolutionMasstransport(IssmDouble* solution);
@@ -257,4 +260,5 @@
 		Seg*	         SpawnSeg(int index1,int index2);
 		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
+		void           TransformSolutionCoord(IssmDouble* values,int transformenum);
 		#ifdef _HAVE_STRESSBALANCE_
 		ElementMatrix* CreateKMatrixStressbalanceSSA(void);