Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16783)
@@ -2274,270 +2274,4 @@
 		this->inputs->AddInput(datasetinput);
 	}
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolution {{{*/
-void  Penta::InputUpdateFromSolution(IssmDouble* solution){
-
-	int analysis_type,inputenum;
-
-	/*retreive parameters: */
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-	switch(analysis_type){
-	#ifdef _HAVE_STRESSBALANCE_
-	case StressbalanceAnalysisEnum:
-		InputUpdateFromSolutionStressbalanceHoriz( solution);
-		break;
-	case StressbalanceSIAAnalysisEnum:
-		InputUpdateFromSolutionStressbalanceSIA( solution);
-		break;
-	case StressbalanceVerticalAnalysisEnum:
-		InputUpdateFromSolutionStressbalanceVert( solution);
-		break;
-	#endif
-	#ifdef _HAVE_CONTROL_
-	case AdjointHorizAnalysisEnum:
-		int approximation;
-		inputs->GetInputValue(&approximation,ApproximationEnum);
-		if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
-			InputUpdateFromSolutionAdjointFS( solution);
-		}
-		else{
-			InputUpdateFromSolutionAdjointHoriz( solution);
-		}
-		break;
-	#endif
-	#ifdef _HAVE_THERMAL_
-	case ThermalAnalysisEnum:
-		InputUpdateFromSolutionThermal( solution);
-		break;
-	case EnthalpyAnalysisEnum:
-		InputUpdateFromSolutionEnthalpy( solution);
-		break;
-	case MeltingAnalysisEnum:
-		InputUpdateFromSolutionOneDof(solution,BasalforcingsMeltingRateEnum);
-		break;
-	#endif
-	case L2ProjectionBaseAnalysisEnum:
-		this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
-		InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
-		break;
-	case MasstransportAnalysisEnum:
-		InputUpdateFromSolutionMasstransport(solution);
-		break;
-	case FreeSurfaceTopAnalysisEnum:
-		InputUpdateFromSolutionFreeSurfaceTop(solution);
-		break;
-	case FreeSurfaceBaseAnalysisEnum:
-		InputUpdateFromSolutionFreeSurfaceBase(solution);
-		break;
-	#ifdef _HAVE_BALANCED_
-	case BalancethicknessAnalysisEnum:
-		InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum);
-		break;
-	#endif
-	#ifdef _HAVE_HYDROLOGY_
-	case HydrologyDCInefficientAnalysisEnum:
-		InputUpdateFromSolutionHydrologyDCInefficient(solution);
-		break;
-	case HydrologyDCEfficientAnalysisEnum:
-		InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
-		break;
-	case L2ProjectionEPLAnalysisEnum:
-		this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
-		InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
-		break;
-	#endif
-	default: 
-		_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionMasstransport{{{*/
-void  Penta::InputUpdateFromSolutionMasstransport(IssmDouble* solution){
-
-	const int  numdof   = NDOF1*NUMVERTICES;
-	const int  numdof2d = NDOF1*NUMVERTICES2D;
-
-	int    i,hydroadjustment;
-	int*   doflist = NULL;
-	IssmDouble rho_ice,rho_water,minthickness;
-	IssmDouble newthickness[numdof];
-	IssmDouble newbed[numdof];
-	IssmDouble newsurface[numdof];
-	IssmDouble oldbed[NUMVERTICES];
-	IssmDouble oldsurface[NUMVERTICES];
-	IssmDouble oldthickness[NUMVERTICES];
-	IssmDouble phi[NUMVERTICES];
-	Penta  *penta   = NULL;
-
-	/*If not on bed, return*/
-	if (!IsOnBed()) return;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector and extrude it */
-	this->parameters->FindParam(&minthickness,MasstransportMinThicknessEnum);
-	for(i=0;i<numdof2d;i++){
-		newthickness[i]=solution[doflist[i]];
-		if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
-		/*Constrain thickness to be at least 1m*/
-		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
-		newthickness[i+numdof2d]=newthickness[i];
-	}
-
-	/*Get previous bed, thickness and surface*/
-	GetInputListOnVertices(&oldbed[0],BedEnum);
-	GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
-	GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
-	GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
-
-	/*Fing MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
-	this->parameters->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
-
-	/*recover material parameters: */
-	rho_ice=matpar->GetRhoIce();
-	rho_water=matpar->GetRhoWater();
-
-	for(i=0;i<numdof;i++) {
-		/*If shelf: hydrostatic equilibrium*/
-		if (phi[i]>0.){
-			newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
-			newbed[i]=oldbed[i];               //same bed: do nothing
-		}
-		else{ //so it is an ice shelf
-			if(hydroadjustment==AbsoluteEnum){
-				newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
-				newbed[i]=newthickness[i]*(-rho_ice/rho_water);
-			}
-			else if(hydroadjustment==IncrementalEnum){
-				newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH 
-				newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
-			}
-			else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
-		}
-	}
-
-	/*Start looping over all elements above current element and update all inputs*/
-	penta=this;
-	for(;;){
-		/*Add input to the element: */
-		penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
-		penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
-		penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if (penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceTop{{{*/
-void  Penta::InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solution){
-
-	const int  numdof   = NDOF1*NUMVERTICES;
-	const int  numdof2d = NDOF1*NUMVERTICES2D;
-
-	int    i;
-	int*   doflist = NULL;
-	IssmDouble newthickness[numdof];
-	IssmDouble newbed[numdof];
-	IssmDouble newsurface[numdof];
-
-	/*If not on bed, return*/
-	if (!IsOnSurface()) return;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector and extrude it */
-	for(i=0;i<numdof2d;i++){
-		newsurface[i+numdof2d]=solution[doflist[i+numdof2d]];
-		if(xIsNan<IssmDouble>(newsurface[i+numdof2d])) _error_("NaN found in solution vector");
-		newsurface[i]=newsurface[i+numdof2d];
-	}
-
-	/*Get previous bed and thickness*/
-	GetInputListOnVertices(&newbed[0],BedEnum);
-
-	for(i=0;i<numdof;i++) {
-		newthickness[i]=newsurface[i]-newbed[i];
-		_assert_(newthickness[i]>0.);
-	}
-
-	/*Start looping over all elements above current element and update all inputs*/
-	Penta* penta=this;
-	for(;;){
-		/*Add input to the element: */
-		penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
-		penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if(penta->IsOnBed()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetLowerElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceBase{{{*/
-void  Penta::InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solution){
-
-	const int  numdof   = NDOF1*NUMVERTICES;
-	const int  numdof2d = NDOF1*NUMVERTICES2D;
-
-	int    i;
-	int*   doflist = NULL;
-	IssmDouble newthickness[numdof];
-	IssmDouble newbed[numdof];
-	IssmDouble newsurface[numdof];
-
-	/*If not on bed, return*/
-	if (!IsOnBed()) return;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector and extrude it */
-	for(i=0;i<numdof2d;i++){
-		newbed[i]=solution[doflist[i]];
-		if(xIsNan<IssmDouble>(newbed[i])) _error_("NaN found in solution vector");
-		newbed[i+numdof2d]=newbed[i];
-	}
-
-	/*Get previous bed and thickness*/
-	GetInputListOnVertices(&newsurface[0],SurfaceEnum);
-
-	for(i=0;i<numdof;i++) {
-		newthickness[i]=newsurface[i]-newbed[i];
-		_assert_(newthickness[i]>0.);
-	}
-
-	/*Start looping over all elements above current element and update all inputs*/
-	Penta* penta=this;
-	for(;;){
-		/*Add input to the element: */
-		penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
-		penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if(penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
 }
 /*}}}*/
@@ -4854,156 +4588,4 @@
 	delete friction;
 	return pe;
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionThermal {{{*/
-void  Penta::InputUpdateFromSolutionThermal(IssmDouble* solution){
-
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	bool        converged;
-	int         i,rheology_law;
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  values[numdof];
-	IssmDouble  B[numdof],surface[numdof];
-	IssmDouble  B_average;
-	int        *doflist = NULL;
-	bool        hack    = false;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++){
-		values[i]=solution[doflist[i]];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-		//if(values[i]<0)      _printf_("temperature < 0°K found in solution vector\n");
-		//if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
-	}
-
-	if(hack){
-		/*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
-		IssmDouble pressure[numdof];
-		GetInputListOnVertices(&pressure[0],PressureEnum);
-		for(i=0;i<numdof;i++){
-			if(values[i]>matpar->TMeltingPoint(pressure[i]))     values[i]=matpar->TMeltingPoint(pressure[i]);
-			if(values[i]<matpar->TMeltingPoint(pressure[i])-50.) values[i]=matpar->TMeltingPoint(pressure[i])-50.;
-		}
-	}
-
-
-	/*Get all inputs and parameters*/
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-
-	this->inputs->GetInputValue(&converged,ConvergedEnum);
-	if(converged){
-		this->inputs->AddInput(new PentaInput(TemperatureEnum,values,P1Enum));
-
-		/*Update Rheology only if converged (we must make sure that the temperature is below melting point
-		 * otherwise the rheology could be negative*/
-		this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
-		GetInputListOnVertices(&surface[0],SurfaceEnum);
-		switch(rheology_law){
-			case NoneEnum:
-				/*Do nothing: B is not temperature dependent*/
-				break;
-			case PatersonEnum:
-				for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
-				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
-				break;
-			case ArrheniusEnum:
-				for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
-				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
-				break;
-			default:
-				_error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
-
-		}
-	}
-	else{
-		this->inputs->AddInput(new PentaInput(TemperaturePicardEnum,values,P1Enum));
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{*/
-void  Penta::InputUpdateFromSolutionEnthalpy(IssmDouble* solution){
-
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	bool       converged=false;
-	int        i,rheology_law;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble values[numdof];
-	IssmDouble pressure[numdof];
-	IssmDouble surface[numdof];
-	IssmDouble temperatures[numdof];
-	IssmDouble waterfraction[numdof];
-	IssmDouble B[numdof];
-	int*       doflist=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++){
-		values[i]=solution[doflist[i]];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Get all inputs and parameters*/
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	GetInputListOnVertices(&pressure[0],PressureEnum);
-	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-
-	this->inputs->GetInputValue(&converged,ConvergedEnum);
-	if(converged){
-		/*Convert enthalpy into temperature and water fraction*/
-		for(i=0;i<numdof;i++){
-			matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
-			if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
-			//if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
-		}
-
-		this->inputs->AddInput(new PentaInput(EnthalpyEnum,values,P1Enum));
-		this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
-		this->inputs->AddInput(new PentaInput(TemperatureEnum,temperatures,P1Enum));
-
-		/*Update Rheology only if converged (we must make sure that the temperature is below melting point
-		 * otherwise the rheology could be negative*/
-		this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
-		GetInputListOnVertices(&surface[0],SurfaceEnum);
-		switch(rheology_law){
-			case NoneEnum:
-				/*Do nothing: B is not temperature dependent*/
-				break;
-			case PatersonEnum:
-				for(i=0;i<numdof;i++) B[i]=Paterson(temperatures[i]);
-				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
-				break;
-			case ArrheniusEnum:
-				for(i=0;i<numdof;i++) B[i]=Arrhenius(temperatures[i],surface[i]-xyz_list[i][2],material->GetN());
-				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
-				break;
-			case LliboutryDuvalEnum:
-				for(i=0;i<numdof;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat());
-				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
-				break;
-			default:
-				_error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
-		}
-	}
-	else{
-		this->inputs->AddInput(new PentaInput(EnthalpyPicardEnum,values,P1Enum));
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
 }
 /*}}}*/
@@ -6371,103 +5953,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionAdjointFS {{{*/
-void  Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){
-
-	int          i;
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	IssmDouble   FSreconditioning;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = this->NumberofNodesVelocity();
-	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
-	int pnumdof   = pnumnodes*NDOF1;
-
-	/*Initialize values*/
-	IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
-	IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
-	IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);
-
-	/*Get dof list: */
-	GetDofListVelocity(&vdoflist,GsetEnum);
-	GetDofListPressure(&pdoflist,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
-	for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
-
-	/*Transform solution in Cartesian Space*/
-	::TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
-
-	/*fill in all arrays: */
-	for(i=0;i<vnumnodes;i++){
-		lambdax[i] = vvalues[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
-		lambday[i] = vvalues[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
-		lambdaz[i] = vvalues[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
-	}
-	for(i=0;i<pnumnodes;i++){
-		lambdap[i] = pvalues[i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Recondition pressure and compute vel: */
-	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
-	this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
-	this->inputs->AddInput(new PentaInput(AdjointzEnum,lambdaz,P1Enum));
-	this->inputs->AddInput(new PentaInput(AdjointpEnum,lambdap,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(vdoflist);
-	xDelete<int>(pdoflist);
-	xDelete<IssmDouble>(lambdap);
-	xDelete<IssmDouble>(lambdaz);
-	xDelete<IssmDouble>(lambday);
-	xDelete<IssmDouble>(lambdax);
-	xDelete<IssmDouble>(vvalues);
-	xDelete<IssmDouble>(pvalues);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionAdjointHoriz {{{*/
-void  Penta::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){
-
-	const int numdof=NDOF2*NUMVERTICES;
-
-	int    i;
-	IssmDouble values[numdof];
-	IssmDouble lambdax[NUMVERTICES];
-	IssmDouble lambday[NUMVERTICES];
-	int*   doflist=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<NUMVERTICES;i++){
-		lambdax[i]=values[i*NDOF2+0];
-		lambday[i]=values[i*NDOF2+1];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(lambdax[i]))       _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(lambday[i]))       _error_("NaN found in solution vector");
-	}
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
-	this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
 /*FUNCTION Penta::SurfaceAverageVelMisfit {{{*/
 IssmDouble Penta::SurfaceAverageVelMisfit(void){
@@ -9846,825 +9329,4 @@
 	*pviscosity = viscosity;
 	return; 
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHoriz {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceHoriz(IssmDouble* solution){
-
-	int  approximation;
-
-	/*Recover inputs*/
-	inputs->GetInputValue(&approximation,ApproximationEnum);
-
-	/*SSA, everything is done by the element on bed*/
-	if (approximation==SSAApproximationEnum){
-		if (!IsOnBed()){
-			/*Do nothing. Element on bed will take care of it*/
-			return;
-		}
-		else{
-			InputUpdateFromSolutionStressbalanceSSA(solution);
-			return;
-		}
-	}
-	if (approximation==L1L2ApproximationEnum){
-		if (!IsOnBed()) return;
-		InputUpdateFromSolutionStressbalanceL1L2(solution);
-		return;
-	}
-	else if (approximation==HOApproximationEnum){
-		InputUpdateFromSolutionStressbalanceHO(solution);
-	}
-	else if (approximation==HOFSApproximationEnum){
-		InputUpdateFromSolutionStressbalanceHOFS(solution);
-	}
-	else if (approximation==SSAFSApproximationEnum){
-		InputUpdateFromSolutionStressbalanceSSAFS(solution);
-	}
-	else if (approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
-		InputUpdateFromSolutionStressbalanceFS(solution);
-	}
-	else if (approximation==SSAHOApproximationEnum){
-		InputUpdateFromSolutionStressbalanceSSAHO(solution);
-	}
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSA {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceSSA(IssmDouble* solution){
-
-	int         numnodes = this->NumberofNodes();
-	int         numdof   = NDOF2*numnodes;
-
-	int         i;
-	IssmDouble  rho_ice,g;
-	IssmDouble  values[2*NUMVERTICES];
-	IssmDouble  vx[NUMVERTICES];
-	IssmDouble  vy[NUMVERTICES];
-	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vel[NUMVERTICES];
-	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  surface[NUMVERTICES];
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int        *doflist = NULL;
-	Penta      *penta   = NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
-
-	/*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*/
-	::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 */
-	for(i=0;i<3;i++){
-		vx[i]  =values[i*NDOF2+0];
-		vy[i]  =values[i*NDOF2+1];
-		vx[i+3]=vx[i];
-		vy[i+3]=vy[i];
-
-		/*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 parameters fro pressure computation*/
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-
-	/*Start looping over all elements above current element and update all inputs*/
-	penta=this;
-	for(;;){
-
-		/*Get node data: */
-		::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
-
-		/*Now Compute vel*/
-		GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
-		for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
-
-		/*Now compute pressure*/
-		GetInputListOnVertices(&surface[0],SurfaceEnum);
-		for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-
-		/*Now, we have to move the previous Vx and Vy inputs  to old 
-		 * status, otherwise, we'll wipe them off: */
-		penta->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-		penta->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-		penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-		/*Add vx and vy as inputs to the tria element: */
-		penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-		penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-		penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-		penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if (penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSAHO {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceSSAHO(IssmDouble* solution){
-
-	const int    numdof=NDOF2*NUMVERTICES;
-	const int    numdof2d=NDOF2*NUMVERTICES2D;
-
-	int     i;
-	IssmDouble  rho_ice,g;
-	IssmDouble  SSA_values[numdof];
-	IssmDouble  HO_values[numdof];
-	IssmDouble  vx[NUMVERTICES];
-	IssmDouble  vy[NUMVERTICES];
-	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vel[NUMVERTICES];
-	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  surface[NUMVERTICES];
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int*    doflistp = NULL;
-	int*    doflistm = NULL;
-	Penta   *penta   = NULL;
-
-	/*OK, we have to add results of this element for HO 
-	 * and results from the penta at base for SSA. Now recover results*/
-	penta=(Penta*)GetBasalElement();
-
-	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
-	GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
-	penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
-
-	/*Get node data: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof2d;i++){
-		HO_values[i]=solution[doflistp[i]];
-		SSA_values[i]=solution[doflistm[i]];
-	}
-	for(i=numdof2d;i<numdof;i++){
-		HO_values[i]=solution[doflistp[i]];
-		SSA_values[i]=SSA_values[i-numdof2d];
-	}
-
-	/*Transform solution in Cartesian Space*/
-	::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: */
-	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=SSA_values[i*NDOF2+0]+HO_values[i*NDOF2+0];
-		vy[i]=SSA_values[i*NDOF2+1]+HO_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");
-	}
-
-	/*Now Compute vel*/
-	GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
-
-	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
-	 *so the pressure is just the pressure at the z elevation: */
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-	GetInputListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflistp);
-	xDelete<int>(doflistm);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSAFS {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceSSAFS(IssmDouble* solution){
-
-	const int    numdofSSA = NDOF2*NUMVERTICES;
-	const int    numdof2d  = NDOF2*NUMVERTICES2D;
-	const int    numdofFSv = NDOF3*NUMVERTICES;
-	const int    numdofFSp = NDOF1*NUMVERTICES;
-
-	int     i;
-	IssmDouble  FSreconditioning;
-	IssmDouble  SSA_values[numdofSSA];
-	IssmDouble  FS_values[numdofFSv+numdofFSp];
-	IssmDouble  vx[NUMVERTICES];
-	IssmDouble  vy[NUMVERTICES];
-	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vzSSA[NUMVERTICES];
-	IssmDouble  vzFS[NUMVERTICES];
-	IssmDouble  vel[NUMVERTICES];
-	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int   *doflistSSA = NULL;
-	int   *doflistFSv = NULL;
-	int   *doflistFSp = NULL;
-	Penta *penta      = NULL;
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(2*NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
-	for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
-
-	/*OK, we have to add results of this element for SSA 
-	 * and results from the penta at base for SSA. Now recover results*/
-	penta=(Penta*)GetBasalElement();
-
-	/*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
-	penta->GetDofList(&doflistSSA,SSAApproximationEnum,GsetEnum);
-	GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
-	GetDofListPressure(&doflistFSp,GsetEnum);
-	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-
-	/*Get node data: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof2d;i++){
-		SSA_values[i]=solution[doflistSSA[i]];
-		SSA_values[i+numdof2d]=solution[doflistSSA[i]];
-	}
-	for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
-	for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
-
-	/*Transform solution in Cartesian Space*/
-	::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: */
-	for(i=0;i<NUMVERTICES;i++){
-		vx[i]       = FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
-		vy[i]       = FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
-		vzFS[i]     = FS_values[i*NDOF3+2];
-		pressure[i] = FS_values[NUMVERTICES*NDOF3+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*/
-	Input* vzSSA_input=inputs->GetInput(VzSSAEnum);
-	if (vzSSA_input){
-		if (vzSSA_input->ObjectEnum()!=PentaInputEnum){
-			_error_("Cannot compute Vel as VzSSA is of type " << EnumToStringx(vzSSA_input->ObjectEnum()));
-		}
-		GetInputListOnVertices(&vzSSA[0],VzSSAEnum);
-	}
-	else{
-		_error_("Cannot update solution as VzSSA is not present");
-	}
-
-	/*Now Compute vel*/
-	for(i=0;i<NUMVERTICES;i++) {
-		vz[i]  = vzSSA[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: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflistSSA);
-	xDelete<int>(doflistFSv);
-	xDelete<int>(doflistFSp);
-	xDelete<int>(cs_list);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceL1L2 {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceL1L2(IssmDouble* solution){
-
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int     i;
-	IssmDouble  rho_ice,g;
-	IssmDouble  values[numdof];
-	IssmDouble  vx[NUMVERTICES];
-	IssmDouble  vy[NUMVERTICES];
-	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vel[NUMVERTICES];
-	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  surface[NUMVERTICES];
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int    *doflist = NULL;
-	Penta  *penta   = NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,L1L2ApproximationEnum,GsetEnum);
-
-	/*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*/
-	::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 */
-	for(i=0;i<3;i++){
-		vx[i]  =values[i*NDOF2+0];
-		vy[i]  =values[i*NDOF2+1];
-		vx[i+3]=vx[i];
-		vy[i+3]=vy[i];
-
-		/*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 parameters fro pressure computation*/
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-
-	/*Start looping over all elements above current element and update all inputs*/
-	penta=this;
-	for(;;){
-
-		/*Get node data: */
-		::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
-
-		/*Now Compute vel*/
-		GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
-		for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
-
-		/*Now compute pressure*/
-		GetInputListOnVertices(&surface[0],SurfaceEnum);
-		for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-
-		/*Now, we have to move the previous Vx and Vy inputs  to old 
-		 * status, otherwise, we'll wipe them off: */
-		penta->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-		penta->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-		penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-		/*Add vx and vy as inputs to the tria element: */
-		penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-		penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-		penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-		penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if (penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHO {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceHO(IssmDouble* solution){
-
-	int         i;
-	IssmDouble  rho_ice,g;
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int        *doflist = NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	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*/
-	::TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*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*/
-	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 3D, 
-	 *so the pressure is just the pressure at the z elevation: */
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-	GetInputListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(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<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHOFS {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceHOFS(IssmDouble* solution){
-
-	const int    numdofHO  = NDOF2*NUMVERTICES;
-	const int    numdofFSv = NDOF3*NUMVERTICES;
-	const int    numdofFSp = NDOF1*NUMVERTICES;
-
-	int        i;
-	IssmDouble HO_values[numdofHO];
-	IssmDouble FS_values[numdofFSv+numdofFSp];
-	IssmDouble vx[NUMVERTICES];
-	IssmDouble vy[NUMVERTICES];
-	IssmDouble vz[NUMVERTICES];
-	IssmDouble vzHO[NUMVERTICES];
-	IssmDouble vzFS[NUMVERTICES];
-	IssmDouble vel[NUMVERTICES];
-	IssmDouble pressure[NUMVERTICES];
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble FSreconditioning;
-	int*       doflistHO        = NULL;
-	int*       doflistFSv        = NULL;
-	int*       doflistFSp = NULL;
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(2*NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
-	for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
-
-	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
-	GetDofList(&doflistHO,HOApproximationEnum,GsetEnum);
-	GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
-	GetDofListPressure(&doflistFSp,GsetEnum);
-	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-
-	/*Get node data: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdofHO;i++)  HO_values[i]=solution[doflistHO[i]];
-	for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
-	for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
-
-	/*Transform solution in Cartesian Space*/
-	::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: */
-	for(i=0;i<NUMVERTICES;i++){
-		vx[i]       = FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
-		vy[i]       = FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
-		vzFS[i]     = FS_values[i*NDOF3+2];
-		pressure[i] = FS_values[NUMVERTICES*NDOF3+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*/
-	Input* vzHO_input=inputs->GetInput(VzHOEnum);
-	if (vzHO_input){
-		if (vzHO_input->ObjectEnum()!=PentaInputEnum){
-			_error_("Cannot compute Vel as VzHO is of type " << EnumToStringx(vzHO_input->ObjectEnum()));
-		}
-		GetInputListOnVertices(&vzHO[0],VzHOEnum);
-	}
-	else{
-		_error_("Cannot update solution as VzHO is not present");
-	}
-
-	/*Now Compute vel*/
-	for(i=0;i<NUMVERTICES;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: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflistHO);
-	xDelete<int>(doflistFSv);
-	xDelete<int>(doflistFSp);
-	xDelete<int>(cs_list);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSIA {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
-
-	int         numnodes = this->NumberofNodes();
-	int         numdof=NDOF2*numnodes;
-
-	int     i;
-	IssmDouble  rho_ice,g;
-	IssmDouble  values[numdof];
-	IssmDouble  vx[NUMVERTICES];
-	IssmDouble  vy[NUMVERTICES];
-	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vel[NUMVERTICES];
-	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  surface[NUMVERTICES];
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	int*    doflist = NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Get node data: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<NUMVERTICES;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");
-	}
-
-	/*Now Compute vel*/
-	GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
-
-	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
-	 *so the pressure is just the pressure at the z elevation: */
-	rho_ice=matpar->GetRhoIce();
-	g=matpar->GetG();
-	GetInputListOnVertices(&surface[0],SurfaceEnum);
-	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceVert {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceVert(IssmDouble* solution){
-
-	int          numnodes = this->NumberofNodes();
-	int          numdof=NDOF1*numnodes;
-
-	int          i;
-	int          approximation;
-	IssmDouble   rho_ice,g;
-	IssmDouble   values[numdof];
-	IssmDouble   vx[NUMVERTICES];
-	IssmDouble   vy[NUMVERTICES];
-	IssmDouble   vz[NUMVERTICES];
-	IssmDouble   vzSSA[NUMVERTICES];
-	IssmDouble   vzHO[NUMVERTICES];
-	IssmDouble   vzFS[NUMVERTICES];
-	IssmDouble   vel[NUMVERTICES];
-	IssmDouble   pressure[NUMVERTICES];
-	IssmDouble   surface[NUMVERTICES];
-	IssmDouble   xyz_list[NUMVERTICES][3];
-	int*         doflist      = NULL;
-
-	/*Get the approximation and do nothing if the element in FS or None*/
-	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
-		return;
-	}
-
-	/*Get dof list and vertices coordinates: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Use the dof list to index into the solution vector vz: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-	for(i=0;i<NUMVERTICES;i++){
-		vz[i]=values[i*NDOF1+0];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Get Vx and Vy*/
-	GetInputListOnVertices(&vx[0],VxEnum,0.0); //default is 0
-	GetInputListOnVertices(&vy[0],VyEnum,0.0); //default is 0
-
-	/*Do some modifications if we actually have a HOFS or SSAFS element*/
-	if(approximation==HOFSApproximationEnum){
-		Input* vzFS_input=inputs->GetInput(VzFSEnum);
-		if (vzFS_input){
-			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
-			GetInputListOnVertices(&vzFS[0],VzFSEnum);
-		}
-		else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
-		for(i=0;i<NUMVERTICES;i++){
-			vzHO[i]=vz[i];
-			vz[i]=vzHO[i]+vzFS[i];
-		}
-	}
-	else if(approximation==SSAFSApproximationEnum){
-		Input* vzFS_input=inputs->GetInput(VzFSEnum);
-		if (vzFS_input){
-			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
-			GetInputListOnVertices(&vzFS[0],VzFSEnum);
-		}
-		else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
-		for(i=0;i<NUMVERTICES;i++){
-			vzSSA[i]=vz[i];
-			vz[i]=vzSSA[i]+vzFS[i];
-		}
-	}
-
-	/*Now Compute vel*/
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
-
-	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
-	 *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
-	if(approximation!=HOFSApproximationEnum &&  approximation!=SSAFSApproximationEnum){
-		rho_ice=matpar->GetRhoIce();
-		g=matpar->GetG();
-		GetInputListOnVertices(&surface[0],SurfaceEnum);
-		for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-	}
-
-	/*Now, we have to move the previous Vz inputs to old 
-	 * status, otherwise, we'll wipe them off and add the new inputs: */
-	this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
-
-	if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
-		this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-		this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
-	}
-	else if(approximation==HOFSApproximationEnum){
-		this->inputs->AddInput(new PentaInput(VzHOEnum,vzHO,P1Enum));
-	}
-	else if(approximation==SSAFSApproximationEnum){
-		this->inputs->AddInput(new PentaInput(VzSSAEnum,vzSSA,P1Enum));
-	}
-	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionStressbalanceFS {{{*/
-void  Penta::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
-
-	int          i;
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	IssmDouble   FSreconditioning;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = this->NumberofNodesVelocity();
-	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF3;
-	int pnumdof   = pnumnodes*NDOF1;
-
-	/*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);
-	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
-	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
-
-	/*Get dof list: */
-	GetDofListVelocity(&vdoflist,GsetEnum);
-	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*/
-	::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
-
-	/*Ok, we have vx and vy in values, fill in all arrays: */
-	for(i=0;i<vnumnodes;i++){
-		vx[i] = values[i*NDOF3+0];
-		vy[i] = values[i*NDOF3+1];
-		vz[i] = values[i*NDOF3+2];
-		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>(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: */
-	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
-	for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-
-	/*Now, we have to move the previous inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,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);
 }
 /*}}}*/
@@ -10952,73 +9614,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
-void  Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
-
-	const int   numdof   = NDOF1*NUMVERTICES;
-	const int   numdof2d = NDOF1*NUMVERTICES2D;
-	int*        doflist  = NULL;
-	bool        converged;
-	IssmDouble  values[numdof];
-	IssmDouble  residual[numdof];
-	IssmDouble  penalty_factor;
-	IssmDouble  kmax, kappa, h_max;
-	Penta      *penta    = NULL;
-
-	/*If not on bed, return*/
-	if (!IsOnBed()) return;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector and extrude it */
-	for(int i=0;i<numdof2d;i++){
-		values[i]         =solution[doflist[i]];
-		values[i+numdof2d]=values[i];
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-	}
-
-	/*If converged keep the residual in mind*/
-	this->inputs->GetInputValue(&converged,ConvergedEnum);
-
-	/*Get inputs*/
-	if(converged){
-		this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
-		this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-
-		kappa=kmax*pow(10.,penalty_factor);
-
-		Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.	
-		for(int i=0;i<NUMVERTICES2D;i++){
-			tria->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
-			if(values[i]>h_max){
-				residual[i]=kappa*(values[i]-h_max);
-				residual[i+numdof2d]=residual[i];
-			}
-			else{
-				residual[i]=0.0;
-				residual[i+numdof2d]=residual[i];
-			}
-		}
-		delete tria->material; delete tria;
-	}
-
-	/*Start looping over all elements above current element and update all inputs*/
-	penta=this;
-	for(;;){
-		/*Add input to the element: */
-		penta->inputs->AddInput(new PentaInput(SedimentHeadEnum,values,P1Enum));
-		penta->inputs->AddInput(new PentaInput(SedimentHeadResidualEnum,residual,P1Enum));
-
-		/*Stop if we have reached the surface*/
-		if (penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
 /*FUNCTION Penta::UpdateConstraintsL2ProjectionEPL{{{*/
 void  Penta::UpdateConstraintsL2ProjectionEPL(void){
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16783)
@@ -60,5 +60,4 @@
 		void  InputUpdateFromConstant(IssmDouble constant, int name);
 		void  InputUpdateFromConstant(int constant, int name);
-		void  InputUpdateFromSolution(IssmDouble* solutiong);
 		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
 		#ifdef _HAVE_DAKOTA_
@@ -232,7 +231,4 @@
 		void           InputChangeName(int input_enum, int enum_type_old);
 		void	         InputExtrude(int enum_type,int object_type);
-		void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
-		void           InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solutiong);
-		void           InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solutiong);
 		void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
 		void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
@@ -290,14 +286,4 @@
 		ElementMatrix* CreateJacobianStressbalanceHO(void);
 		ElementMatrix* CreateJacobianStressbalanceFS(void);
-		void           InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceSSA( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceSSAHO( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceSSAFS( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceL1L2( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceHO( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceHOFS( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceVert( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionStressbalanceFS( IssmDouble* solutiong);
 		ElementVector* CreatePVectorCouplingSSAFS(void);
 		ElementVector* CreatePVectorCouplingSSAFSViscous(void);
@@ -335,6 +321,4 @@
 		ElementVector* CreatePVectorAdjointHO(void);
 		ElementVector* CreatePVectorAdjointFS(void);
-		void           InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionAdjointFS( IssmDouble* solutiong);
 		#endif
 
@@ -351,5 +335,4 @@
 		void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
 		void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
-		void           InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
 		void           ComputeEPLThickness(void);
 		void           UpdateConstraintsL2ProjectionEPL(void);
@@ -375,6 +358,4 @@
 		ElementVector* CreatePVectorThermalShelf(void);
 		ElementVector* CreatePVectorThermalSheet(void);
-		void           InputUpdateFromSolutionThermal(IssmDouble* solutiong);
-		void           InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong);
 		void           UpdateBasalConstraintsEnthalpy(void);
 		void           ComputeBasalMeltingrate(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16783)
@@ -54,5 +54,4 @@
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void  InputUpdateFromSolution(IssmDouble* solution){_error_("not implemented yet");};
 		void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
 		void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16783)
@@ -1809,95 +1809,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolution {{{*/
-void  Tria::InputUpdateFromSolution(IssmDouble* solution){
-
-	/*retrive parameters: */
-	int analysis_type,extrusioninput;
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-	switch(analysis_type){
-		#ifdef _HAVE_STRESSBALANCE_
-		case StressbalanceAnalysisEnum:
-			int approximation;
-			inputs->GetInputValue(&approximation,ApproximationEnum);
-			if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
-				InputUpdateFromSolutionStressbalanceFS(solution);
-			}
-			else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){
-				InputUpdateFromSolutionStressbalanceHoriz(solution);
-			}
-			else{
-				_error_("approximation not supported yet");
-			}
-			break;
-		case StressbalanceSIAAnalysisEnum:
-			InputUpdateFromSolutionStressbalanceHoriz(solution);
-			break;
-		#endif
-		#ifdef _HAVE_CONTROL_
-		case AdjointHorizAnalysisEnum:
-			InputUpdateFromSolutionAdjointHoriz(solution);
-			break;
-		case AdjointBalancethicknessAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,AdjointEnum);
-			break;
-		#endif
-		#ifdef _HAVE_HYDROLOGY_ 
-		case HydrologyShreveAnalysisEnum:
-			InputUpdateFromSolutionHydrologyShreve(solution);
-			break;
-		case HydrologyDCInefficientAnalysisEnum:
-			InputUpdateFromSolutionHydrologyDCInefficient(solution);
-			break;
-		case HydrologyDCEfficientAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
-			break;
-		case L2ProjectionEPLAnalysisEnum:
-			this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum);
-			InputUpdateFromSolutionOneDof(solution,extrusioninput);
-			break;
-		#endif
-	 	#ifdef _HAVE_DAMAGE_
-		case DamageEvolutionAnalysisEnum:
-			InputUpdateFromSolutionDamageEvolution(solution);
-			break;
-		#endif
-	 	#ifdef _HAVE_BALANCED_
-		case BalancethicknessAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
-			break;
-		case BalancevelocityAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,VelEnum);
-			break;
-		case SmoothedSurfaceSlopeXAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
-			break;
-		case SmoothedSurfaceSlopeYAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
-			break;
-		#endif
-		case L2ProjectionBaseAnalysisEnum:
-			this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum);
-			InputUpdateFromSolutionOneDof(solution,extrusioninput);
-			break;
-		case MasstransportAnalysisEnum:
-			InputUpdateFromSolutionMasstransport(solution);
-			break;
-		case FreeSurfaceTopAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
-			break;
-		case FreeSurfaceBaseAnalysisEnum:
-			InputUpdateFromSolutionOneDof(solution,BedEnum);
-			break;
-		case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum:
-			this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
-			InputUpdateFromSolutionOneDof(solution,extrusioninput);
-			break;
-		default:
-			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
-	}
-}
-/*}}}*/
 /*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{*/
 void  Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
@@ -1924,80 +1833,4 @@
 	/*Free ressources:*/
 	xDelete<IssmDouble>(values);
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionMasstransport{{{*/
-void  Tria::InputUpdateFromSolutionMasstransport(IssmDouble* solution){
-
-	/*Intermediaries*/
-	int        i,hydroadjustment;
-	int*       doflist=NULL;
-	IssmDouble rho_ice,rho_water,minthickness;
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Fetch dof list and allocate solution vector*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
-	IssmDouble* newbed       = xNew<IssmDouble>(numnodes);
-	IssmDouble* newsurface   = xNew<IssmDouble>(numnodes);
-	IssmDouble* oldthickness = xNew<IssmDouble>(numnodes);
-	IssmDouble* oldbed       = xNew<IssmDouble>(numnodes);
-	IssmDouble* oldsurface   = xNew<IssmDouble>(numnodes);
-	IssmDouble* phi          = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	this->parameters->FindParam(&minthickness,MasstransportMinThicknessEnum);
-	for(i=0;i<numnodes;i++){
-		newthickness[i]=solution[doflist[i]];
-		if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
-		/*Constrain thickness to be at least 1m*/
-		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
-	}
-
-	/*Get previous bed, thickness and surface*/
-	GetInputListOnNodes(&oldbed[0],BedEnum);
-	GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
-	GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
-	GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
-
-	/*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
-	this->parameters->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
-	rho_ice=matpar->GetRhoIce();
-	rho_water=matpar->GetRhoWater();
-
-	for(i=0;i<numnodes;i++) {
-		/*If shelf: hydrostatic equilibrium*/
-		if (phi[i]>0.){
-			newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
-			newbed[i]     = oldbed[i];                 //same bed: do nothing
-		}
-		else{ //this is an ice shelf
-			if(hydroadjustment==AbsoluteEnum){
-				newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
-				newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
-			}
-			else if(hydroadjustment==IncrementalEnum){
-				newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
-				newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
-			}
-			else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
-		}
-	}
-
-	/*Add input to the element: */
-	this->inputs->AddInput(new TriaInput(ThicknessEnum,newthickness,P1Enum));
-	this->inputs->AddInput(new TriaInput(SurfaceEnum,newsurface,P1Enum));
-	this->inputs->AddInput(new TriaInput(BedEnum,newbed,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(newthickness);
-	xDelete<IssmDouble>(newbed);
-	xDelete<IssmDouble>(newsurface);
-	xDelete<IssmDouble>(oldthickness);
-	xDelete<IssmDouble>(oldbed);
-	xDelete<IssmDouble>(oldsurface);
-	xDelete<IssmDouble>(phi);
 	xDelete<int>(doflist);
 }
@@ -4414,224 +4247,4 @@
 
 	return y;
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceHoriz {{{*/
-void  Tria::InputUpdateFromSolutionStressbalanceHoriz(IssmDouble* solution){
-
-	int        i;
-	IssmDouble rho_ice,g;
-	int*       doflist=NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	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*/
-	::TransformSolutionCoord(&values[0],nodes,numnodes,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*/
-	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=matpar->GetRhoIce();
-	g=matpar->GetG();
-	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: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new TriaInput(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);
-
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/
-void  Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
-
-	int          i;
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	IssmDouble   FSreconditioning;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = this->NumberofNodesVelocity();
-	int pnumnodes = this->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*NDOF2;
-	int pnumdof   = pnumnodes*NDOF1;
-
-	/*Initialize values*/
-	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
-	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
-	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
-
-	/*Get dof list: */
-	GetDofListVelocity(&vdoflist,GsetEnum);
-	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*/
-	::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
-
-	/*Ok, we have vx and vy in values, fill in all arrays: */
-	for(i=0;i<vnumnodes;i++){
-		vx[i] = values[i*NDOF2+0];
-		vy[i] = values[i*NDOF2+1];
-		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vy[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: */
-	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
-	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: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(vel);
-	xDelete<IssmDouble>(vy);
-	xDelete<IssmDouble>(vx);
-	xDelete<IssmDouble>(values);
-	xDelete<int>(vdoflist);
-	xDelete<int>(pdoflist);
-	xDelete<int>(cs_list);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceSIA {{{*/
-void  Tria::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
-
-	int        i;
-	IssmDouble rho_ice,g;
-	int*       doflist=NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numdof);
-	IssmDouble* vx        = xNew<IssmDouble>(numdof);
-	IssmDouble* vy        = xNew<IssmDouble>(numdof);
-	IssmDouble* vz        = xNew<IssmDouble>(numdof);
-	IssmDouble* vel       = xNew<IssmDouble>(numdof);
-	IssmDouble* pressure  = xNew<IssmDouble>(numdof);
-	IssmDouble* thickness = xNew<IssmDouble>(numdof);
-
-	/*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*/
-	::TransformSolutionCoord(&values[0],nodes,numnodes,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*/
-	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=matpar->GetRhoIce();
-	g=matpar->GetG();
-	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: */
-	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
-	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
-	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
-	this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
-	this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
-	this->inputs->AddInput(new TriaInput(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);
 }
 /*}}}*/
@@ -6340,47 +5953,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{*/
-void  Tria::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){
-
-	int  i;
-	int* doflist=NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values  = xNew<IssmDouble>(numdof);
-	IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
-	IssmDouble* lambday = 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*/
-	::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numnodes;i++){
-		lambdax[i]=values[i*NDOF2+0];
-		lambday[i]=values[i*NDOF2+1];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new TriaInput(AdjointxEnum,lambdax,P1Enum));
-	this->inputs->AddInput(new TriaInput(AdjointyEnum,lambday,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(lambdax);
-	xDelete<IssmDouble>(lambday);
-	xDelete<int>(doflist);
-}
-/*}}}*/
 /*FUNCTION Tria::GetVectorFromControlInputs{{{*/
 void  Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){
@@ -7113,86 +6683,4 @@
 	xDelete<IssmDouble>(values);
 	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/
-void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){
-
-	/*Intermediary*/
-	int* doflist = NULL;
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Fetch dof list and allocate solution vector*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numnodes;i++){
-		values[i]=solution[doflist[i]];
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-		if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
-	}
-
-	/*Add input to the element: */
-	this->inputs->AddInput(new TriaInput(WatercolumnEnum,values,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(values);
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
-void  Tria::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
-
-	/*Intermediaries*/
-	int        *doflist  = NULL;
-	bool        converged;
-	IssmDouble  penalty_factor, dt;
-	IssmDouble  kmax, kappa, h_max;
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Fetch dof list and allocate solution vector*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values   = xNew<IssmDouble>(numnodes);
-	IssmDouble* residual = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numnodes;i++){
-		values[i]=solution[doflist[i]];
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-	}
-
-	/*If converged keep the residual in mind*/
-	this->inputs->GetInputValue(&converged,ConvergedEnum);
-
-	/*Get inputs*/
-	if(converged){
-		this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
-		this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-
-		kappa=kmax*pow(10.,penalty_factor);
-
-		for(int i=0;i<numnodes;i++){
-			this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
-			if(values[i]>h_max){
-				residual[i]=kappa*(values[i]-h_max);
-			}
-			else{
-				residual[i]=0.0;
-			}
-		}
-		this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual,P1Enum));
-	}
-
-	/*Add input to the element: */
-	this->inputs->AddInput(new TriaInput(SedimentHeadEnum,values,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(residual);
-	xDelete<int>(doflist);
 }
 /*}}}*/
@@ -8567,37 +8055,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionDamageEvolution {{{*/
-void  Tria::InputUpdateFromSolutionDamageEvolution(IssmDouble* solution){
-
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	int         i;
-	IssmDouble  values[numdof];
-	IssmDouble  max_damage;
-	int			*doflist = NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Get user-supplied max_damage: */
-	this->parameters->FindParam(&max_damage,DamageMaxDamageEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++){
-		values[i]=solution[doflist[i]];
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-		/*Enforce D < max_damage and D > 0 */
-		if(values[i]>max_damage) values[i]=max_damage;
-		else if(values[i]<0) values[i]=0;
-	}
-
-	/*Get all inputs and parameters*/
-	this->material->inputs->AddInput(new TriaInput(DamageDbarEnum,values,P1Enum));
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16783)
@@ -55,5 +55,4 @@
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void  InputUpdateFromSolution(IssmDouble* solutiong);
 		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
 		#ifdef _HAVE_DAKOTA_
@@ -266,5 +265,4 @@
 		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
 		void	         InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
-		void	         InputUpdateFromSolutionMasstransport(IssmDouble* solution);
 		bool	         IsInput(int name);
 		Gauss*         NewGauss(void);
@@ -295,7 +293,4 @@
 		ElementMatrix* CreateJacobianStressbalanceSSA(void);
 		IssmDouble     GetYcoord(GaussTria* gauss);
-		void	         InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
-		void	         InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
-		void	         InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
 		#endif
 
@@ -305,5 +300,4 @@
 		ElementVector* CreatePVectorAdjointHoriz(void);
 		ElementVector* CreatePVectorAdjointBalancethickness(void);
-		void	         InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution);
 		#endif
 
@@ -326,9 +320,4 @@
 		ElementVector* CreatePVectorL2ProjectionEPL(void);
 		void           CreateHydrologyWaterVelocityInput(void);
-		void	         InputUpdateFromSolutionHydrology(IssmDouble* solution);
-		void	         InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
-		void           InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
-		void	         InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
-		void	         InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
 		void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
 		void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
@@ -348,5 +337,4 @@
 		ElementVector* CreatePVectorDamageEvolution_CG(void);
 		void           DamageEvolutionF(IssmDouble* flist);
-		void	         InputUpdateFromSolutionDamageEvolution(IssmDouble* solution);
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 16783)
@@ -54,5 +54,4 @@
 		void InputUpdateFromConstant(int constant, int name){/*Do nothing*/};
 		void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");}
-		void InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
 		void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16783)
@@ -407,9 +407,4 @@
 }
 /*}}}*/
-/*FUNCTION Pengrid::InputUpdateFromSolution{{{*/
-void  Pengrid::InputUpdateFromSolution(IssmDouble* solution){
-	/*Nothing updated yet*/
-}
-/*}}}*/		
 
 /*Pengrid management:*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16783)
@@ -64,5 +64,4 @@
 		void  InputUpdateFromConstant(int constant, int name);
 		void  InputUpdateFromConstant(bool constant, int name);
-		void  InputUpdateFromSolution(IssmDouble* solution);
 		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Penpair.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.h	(revision 16783)
@@ -44,5 +44,4 @@
 		void  InputUpdateFromConstant(int constant, int name);
 		void  InputUpdateFromConstant(bool constant, int name);
-		void  InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
 		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 16783)
@@ -65,6 +65,5 @@
 		void    InputUpdateFromConstant(int constant, int name){_error_("Not implemented yet!");}
 		void    InputUpdateFromConstant(bool constant, int name);
-		void    InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
-		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void    InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
 		/*}}}*/
 		/*Load virtual functions definitions: {{{*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 16783)
@@ -854,9 +854,4 @@
 }
 /*}}}*/
-/*FUNCTION Matice::InputUpdateFromSolution{{{*/
-void  Matice::InputUpdateFromSolution(IssmDouble* solution){
-	/*Nothing updated yet*/
-}
-/*}}}*/
 /*FUNCTION Matice::InputUpdateFromIoModel{{{*/
 void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 16783)
@@ -45,5 +45,4 @@
 		void  InputUpdateFromConstant(int constant, int name);
 		void  InputUpdateFromConstant(bool constant, int name);
-		void  InputUpdateFromSolution(IssmDouble* solution);
 		void  InputUpdateFromIoModel(int index, IoModel* iomodel);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16783)
@@ -220,9 +220,4 @@
 }
 /*}}}*/
-/*FUNCTION Matpar::InputUpdateFromSolution{{{*/
-void   Matpar::InputUpdateFromSolution(IssmDouble* solution){
-	/*Nothing updated yet*/
-}
-/*}}}*/
 
 /*Matpar management: */
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16783)
@@ -77,5 +77,4 @@
 		void   InputUpdateFromConstant(int constant, int name);
 		void   InputUpdateFromConstant(bool constant, int name);
-		void   InputUpdateFromSolution(IssmDouble* solution);
 		void   InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Update.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Update.h	(revision 16782)
+++ /issm/trunk-jpl/src/c/classes/Update.h	(revision 16783)
@@ -22,5 +22,4 @@
 		virtual void  InputUpdateFromConstant(int constant, int name)=0;
 		virtual void  InputUpdateFromConstant(bool constant, int name)=0;
-		virtual void  InputUpdateFromSolution(IssmDouble* solution)=0;
 		virtual void  InputUpdateFromIoModel(int index, IoModel* iomodel)=0;
 
Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp	(revision 16782)
+++ /issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp	(revision 16783)
@@ -28,6 +28,6 @@
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		//analysis->InputUpdateFromSolution(solution,element);
-		element->InputUpdateFromSolution(solution);
+		analysis->InputUpdateFromSolution(solution,element);
+		//element->InputUpdateFromSolution(solution);
 	}
 	delete analysis;