Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16772)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16773)
@@ -943,7 +943,4 @@
 			InputUpdateFromSolutionHO(solution,element);
 			return;
-		case SIAApproximationEnum:
-			InputUpdateFromSolutionHO(solution,element);
-			return;
 		case L1L2ApproximationEnum:
 			InputUpdateFromSolutionSSA(solution,element);
@@ -961,5 +958,175 @@
 	}
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
+
+	int          i,dim,meshtype;
+	int*         vdoflist=NULL;
+	int*         pdoflist=NULL;
+	IssmDouble   FSreconditioning;
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+	int vnumdof   = vnumnodes*dim;
+	int pnumdof   = pnumnodes*1;
+
+	/*Initialize values*/
+	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
+	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	if(dim==2){
+		for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+	}
+	else{
+		for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+	}
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Get dof list: */
+	element->GetDofListVelocity(&vdoflist,GsetEnum);
+	element->GetDofListPressure(&pdoflist,GsetEnum);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
+	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	element->TransformSolutionCoord(&values[0],&cs_list[0]);
+
+	/*Ok, we have vx and vy in values, fill in all arrays: */
+	for(i=0;i<vnumnodes;i++){
+		vx[i] = values[i*dim+0];
+		vy[i] = values[i*dim+1];
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+
+		if(dim==3){
+			vz[i] = values[i*dim+2];
+			if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+		}
+	}
+	for(i=0;i<pnumnodes;i++){
+		pressure[i] = values[vnumdof+i];
+		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Recondition pressure and compute vel: */
+	for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+	if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+	else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
+
+	/*Now, we have to move the previous inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+	element->InputChangeName(PressureEnum,PressurePicardEnum);
+	if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
+
+	/*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);
+	if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<int>(vdoflist);
+	xDelete<int>(pdoflist);
+	xDelete<int>(cs_list);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i,meshtype;
+	IssmDouble  rho_ice,g;
+	int*        doflist=NULL;
+	IssmDouble* xyz_list=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,HOApproximationEnum,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* surface   = 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);
+	element->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*/
+	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>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,meshtype;
@@ -1060,5 +1227,5 @@
 	if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,meshtype;
@@ -1066,11 +1233,50 @@
 	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->GetInputListOnVertices(thickness,ThicknessEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+			break;
+		case Mesh3DEnum:
+			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]);
+			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()){xDelete<IssmDouble>(xyz_list); 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 = element->GetNumberOfNodes();
+	int numnodes = basalelement->GetNumberOfNodes();
 	int numdof   = numnodes*2;
 
 	/*Fetch dof list and allocate solution vectors*/
-	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numdof);
 	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
@@ -1078,6 +1284,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: */
@@ -1085,6 +1289,6 @@
 
 	/*Transform solution in Cartesian Space*/
-	element->TransformSolutionCoord(&values[0],XYEnum);
-	element->FindParam(&meshtype,MeshTypeEnum);
+	basalelement->TransformSolutionCoord(&values[0],XYEnum);
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -1099,14 +1303,6 @@
 
 	/*Get Vz and compute vel*/
-	element->GetInputListOnNodes(&vz[0],VzEnum,0.);
+	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]);
-
-	/*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 
@@ -1114,15 +1310,11 @@
 	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(VxEnum,vx,P1Enum);
+	element->AddBasalInput(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);
@@ -1132,101 +1324,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
-}/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
-
-	int          i,dim,meshtype;
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	IssmDouble   FSreconditioning;
-
-	element->FindParam(&meshtype,MeshTypeEnum);
-	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->GetNumberOfNodesVelocity();
-	int pnumnodes = element->GetNumberOfNodesPressure();
-	int vnumdof   = vnumnodes*dim;
-	int pnumdof   = pnumnodes*1;
-
-	/*Initialize values*/
-	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
-	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-	if(dim==2){
-		for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
-	}
-	else{
-		for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
-	}
-	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
-
-	/*Get dof list: */
-	element->GetDofListVelocity(&vdoflist,GsetEnum);
-	element->GetDofListPressure(&pdoflist,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
-	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
-
-	/*Transform solution in Cartesian Space*/
-	element->TransformSolutionCoord(&values[0],&cs_list[0]);
-
-	/*Ok, we have vx and vy in values, fill in all arrays: */
-	for(i=0;i<vnumnodes;i++){
-		vx[i] = values[i*dim+0];
-		vy[i] = values[i*dim+1];
-		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-
-		if(dim==3){
-			vz[i] = values[i*dim+2];
-			if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
-		}
-	}
-	for(i=0;i<pnumnodes;i++){
-		pressure[i] = values[vnumdof+i];
-		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Recondition pressure and compute vel: */
-	for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
-	if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-	else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
-
-	/*Now, we have to move the previous inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	element->InputChangeName(VxEnum,VxPicardEnum);
-	element->InputChangeName(VyEnum,VyPicardEnum);
-	element->InputChangeName(PressureEnum,PressurePicardEnum);
-	if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
-
-	/*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);
-	if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(vel);
-	xDelete<IssmDouble>(vz);
-	xDelete<IssmDouble>(vy);
-	xDelete<IssmDouble>(vx);
-	xDelete<IssmDouble>(values);
-	xDelete<int>(vdoflist);
-	xDelete<int>(pdoflist);
-	xDelete<int>(cs_list);
+	if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
@@ -1252,5 +1348,5 @@
 	/*Fetch dof list and allocate solution vectors*/
 	basalelement->GetDofList(&SSAdoflist,SSAApproximationEnum,GsetEnum);
-	element->GetDofList(     &HOdoflist, HOApproximationEnum, GsetEnum);
+	element     ->GetDofList(&HOdoflist, HOApproximationEnum, GsetEnum);
 	IssmDouble* HOvalues  = xNew<IssmDouble>(numdof);
 	IssmDouble* SSAvalues = xNew<IssmDouble>(numdof);
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16772)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16773)
@@ -24,4 +24,5 @@
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);