Index: issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 24221)
+++ issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 24335)
@@ -659,5 +659,5 @@
 	return numdofs;
 }/*}}}*/
-void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+void StressbalanceAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	/*Intermediaries*/
@@ -730,5 +730,11 @@
 		if(iomodel->my_elements[i]){
 			Element* element=(Element*)elements->GetObjectByOffset(counter);
-			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
+			element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
+
+			/*Need to know the type of approximation for this element*/
+			if(iomodel->Data("md.flowequation.element_equation")){
+				inputs2->SetInput(ApproximationEnum,counter,IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[i])));
+			}
+
 			counter++;
 		}
@@ -736,29 +742,29 @@
 
 	/*Create inputs: */
-	iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
-	iomodel->FetchDataToInput(elements,"md.geometry.surface",SurfaceEnum);
-	iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
-	iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
-	iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
-	iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
-	iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum,0.);
-	iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum,0.);
-	iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcex",LoadingforceXEnum);
-	iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcey",LoadingforceYEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.slr.sealevel",SealevelEnum,0);
+	iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum,0.);
+	iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum,0.);
+	iomodel->FetchDataToInput(inputs2,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.stressbalance.loadingforcey",LoadingforceYEnum);
 	#ifdef LATERALFRICTION
-	iomodel->FetchDataToInput(elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
+	iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
 	#endif
 
 	if(iomodel->domaintype!=Domain2DhorizontalEnum){
-		iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
-		iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
+		iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
+		iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
 	}
 	if(iomodel->domaintype==Domain3DEnum){
-		iomodel->FetchDataToInput(elements,"md.flowequation.borderFS",FlowequationBorderFSEnum);
-		iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcez",LoadingforceZEnum);
-		iomodel->FetchDataToInput(elements,"md.initialization.vz",VzEnum,0.);
+		iomodel->FetchDataToInput(inputs2,elements,"md.flowequation.borderFS",FlowequationBorderFSEnum);
+		iomodel->FetchDataToInput(inputs2,elements,"md.stressbalance.loadingforcez",LoadingforceZEnum);
+		iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vz",VzEnum,0.);
 	}
 	if(isFS){
-		iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum,0.);
+		iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum,0.);
 
 		/*Add basal forcings to compute melt rate*/
@@ -767,5 +773,5 @@
 		switch(basalforcing_model){
 			case FloatingMeltRateEnum:
-				iomodel->FetchDataToInput(elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
 				break;
 			case LinearFloatingMeltRateEnum:
@@ -776,17 +782,17 @@
 				break;
 			case SpatialLinearFloatingMeltRateEnum:
-				iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
-				iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
-				iomodel->FetchDataToInput(elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
 				break;
 			case BasalforcingsPicoEnum:
-				iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
 				break;
 			case BasalforcingsIsmip6Enum:
-				iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
 				break;
 			case BeckmannGoosseFloatingMeltRateEnum:
-				iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
-				iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
 				break;
 			default:
@@ -797,5 +803,5 @@
 	iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS");
 	if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
-		InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
+		InputUpdateFromConstantx(inputs2,elements,0.,SigmaNNEnum);
 	}
 
@@ -804,83 +810,83 @@
 		case 1:
 			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
-			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
 			if(FrictionCoupling==3){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
 			else if(FrictionCoupling==4){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
 			}
 			break;
 		case 2:
-			iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
 			break;
 		case 3:
 			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
-			iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.As",FrictionAsEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
 			if(FrictionCoupling==3){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
 			else if(FrictionCoupling==4){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
 			}
 			break;
 		case 4:
-			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
-			iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
-			iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
 			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
 			break;
 		case 5:
-			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.water_layer",FrictionWaterLayerEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.water_layer",FrictionWaterLayerEnum);
 			break;
 		case 6:
-			iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
-			iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
-			iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
 			break;
 		case 7:
 			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
-			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
 			if(FrictionCoupling==3){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
 			else if(FrictionCoupling==4){
-				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
+				iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
 
 			}
 			break;
 		case 9:
-			iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
-			InputUpdateFromConstantx(elements,1.,FrictionPEnum);
-			InputUpdateFromConstantx(elements,1.,FrictionQEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
+			InputUpdateFromConstantx(inputs2,elements,1.,FrictionPEnum);
+			InputUpdateFromConstantx(inputs2,elements,1.,FrictionQEnum);
 			break;
 		case 10:
-			iomodel->FetchDataToInput(elements,"md.friction.till_friction_angle",FrictionTillFrictionAngleEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.sediment_compressibility_coefficient",FrictionSedimentCompressibilityCoefficientEnum);
-			iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn_max",HydrologyWatercolumnMaxEnum);
-			iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.till_friction_angle",FrictionTillFrictionAngleEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.sediment_compressibility_coefficient",FrictionSedimentCompressibilityCoefficientEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.watercolumn_max",HydrologyWatercolumnMaxEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
 			break;
 		case 11:
-			iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.Cmax",FrictionCmaxEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.Cmax",FrictionCmaxEnum);
 			break;
 		case 12:
-			iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
-			iomodel->FetchDataToInput(elements,"md.friction.f",FrictionfEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
+			iomodel->FetchDataToInput(inputs2,elements,"md.friction.f",FrictionfEnum);
 			break;
 		default:
@@ -889,5 +895,5 @@
 
 #ifdef _HAVE_ANDROID_
-	elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
+	inputs2->DuplicateInput(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
 #endif
 
@@ -1059,5 +1065,5 @@
 
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case FSApproximationEnum:
@@ -1072,5 +1078,5 @@
 
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case SSAApproximationEnum:
@@ -1088,5 +1094,5 @@
 ElementMatrix* StressbalanceAnalysis::CreateKMatrix(Element* element){/*{{{*/
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case SIAApproximationEnum:
@@ -1115,5 +1121,5 @@
 
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case SIAApproximationEnum:
@@ -1142,5 +1148,5 @@
 
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case FSApproximationEnum: case NoneApproximationEnum:
@@ -1178,5 +1184,5 @@
 	int numnodes = element->GetNumberOfNodes();
 	int numdof   = numnodes*dofpernode;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 
 	/*Fetch dof list and allocate solution vector*/
@@ -1185,7 +1191,7 @@
 
 	/*Get inputs*/
-	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=NULL;
-	if(domaintype!=Domain2DverticalEnum){vy_input=element->GetInput(VyEnum); _assert_(vy_input);}
+	Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
+	Input2* vy_input=NULL;
+	if(domaintype!=Domain2DverticalEnum){vy_input=element->GetInput2(VyEnum); _assert_(vy_input);}
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -1216,5 +1222,5 @@
 
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	switch(approximation){
 		case FSApproximationEnum: case NoneApproximationEnum:
@@ -1287,7 +1293,7 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* thickness_input = basalelement->GetInput(ThicknessEnum);_assert_(thickness_input);
-	Input* vx_input        = basalelement->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input        = basalelement->GetInput(VyEnum);       _assert_(vy_input);
+	Input2* thickness_input = basalelement->GetInput2(ThicknessEnum);_assert_(thickness_input);
+	Input2* vx_input        = basalelement->GetInput2(VxEnum);       _assert_(vx_input);
+	Input2* vy_input        = basalelement->GetInput2(VyEnum);       _assert_(vy_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1406,6 +1412,6 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum);
-	Input* surface_input    = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* gllevelset_input = NULL;
+	Input2* surface_input    = element->GetInput2(SurfaceEnum); _assert_(surface_input);
+	Input2* gllevelset_input = NULL;
 
 	/*build friction object, used later on: */
@@ -1415,5 +1421,5 @@
 	if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list);
 	if(friction_style==SubelementFriction2Enum){
-		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		gllevelset_input=element->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
 		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
 	   gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
@@ -1501,5 +1507,5 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.);
-	Input* icelevelset_input = element->GetInput(MaskIceLevelsetEnum); _assert_(icelevelset_input);
+	Input2* icelevelset_input = element->GetInput2(MaskIceLevelsetEnum); _assert_(icelevelset_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1563,9 +1569,9 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
-	Input* vy_input    = NULL;
+	Input2* thickness_input=element->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2* vx_input=element->GetInput2(VxEnum);               _assert_(vx_input);
+	Input2* vy_input    = NULL;
 	if(dim==2){
-		vy_input    = element->GetInput(VyEnum);       _assert_(vy_input);
+		vy_input    = element->GetInput2(VyEnum);       _assert_(vy_input);
 	}
 
@@ -1680,6 +1686,6 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input2*     thickness_input=element->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2*     surface_input  =element->GetInput2(SurfaceEnum);   _assert_(surface_input);
 	IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum);
 
@@ -1744,7 +1750,7 @@
 
 	/*Retrieve all inputs and parameters*/
-	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
-	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2* base_input       = element->GetInput2(BaseEnum);       _assert_(base_input);
+	Input2* sealevel_input       = element->GetInput2(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
@@ -1941,5 +1947,5 @@
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(PressureEnum,pressure,P1Enum);
 	xDelete<IssmDouble>(pressure);
 	xDelete<IssmDouble>(thickness);
@@ -1964,9 +1970,9 @@
 	/*Fetch dof list and allocate solution vectors*/
 	basalelement->GetDofListLocal(&doflist,SSAApproximationEnum,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* values = xNew<IssmDouble>(numdof);
+	IssmDouble* vx     = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy     = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz     = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel    = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
@@ -2000,7 +2006,7 @@
 
 	/*Add vx and vy as inputs to the tria element: */
-	element->AddBasalInput(VxEnum,vx,element->GetElementType());
-	if(dim==2)element->AddBasalInput(VyEnum,vy,element->GetElementType());
-	element->AddBasalInput(VelEnum,vel,element->GetElementType());
+	element->AddBasalInput2(VxEnum,vx,element->GetElementType());
+	if(dim==2)element->AddBasalInput2(VyEnum,vy,element->GetElementType());
+	element->AddBasalInput2(VelEnum,vel,element->GetElementType());
 
 	/*Free ressources:*/
@@ -2062,7 +2068,7 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* vx_input      = element->GetInput(VxEnum);      _assert_(vx_input);
-	Input* vy_input      = element->GetInput(VyEnum);      _assert_(vy_input);
+	Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
+	Input2* vx_input      = element->GetInput2(VxEnum);      _assert_(vx_input);
+	Input2* vy_input      = element->GetInput2(VyEnum);      _assert_(vy_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2149,6 +2155,6 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input2*     thickness_input=element->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2*     surface_input  =element->GetInput2(SurfaceEnum);   _assert_(surface_input);
 	IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum);
 
@@ -2200,7 +2206,7 @@
 
 	/*Retrieve all inputs and parameters*/
-	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
-	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2* base_input       = element->GetInput2(BaseEnum);       _assert_(base_input);
+	Input2* sealevel_input       = element->GetInput2(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
@@ -2270,5 +2276,5 @@
 		for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
 	}
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(PressureEnum,pressure,P1Enum);
 	xDelete<IssmDouble>(pressure);
 	xDelete<IssmDouble>(thickness);
@@ -2323,7 +2329,7 @@
 
 	/*Add vx and vy as inputs to the tria element: */
-	element->AddBasalInput(VxEnum,vx,element->GetElementType());
-	element->AddBasalInput(VyEnum,vy,element->GetElementType());
-	element->AddBasalInput(VelEnum,vel,element->GetElementType());
+	element->AddBasalInput2(VxEnum,vx,element->GetElementType());
+	element->AddBasalInput2(VyEnum,vy,element->GetElementType());
+	element->AddBasalInput2(VelEnum,vel,element->GetElementType());
 
 	/*Free ressources:*/
@@ -2359,6 +2365,6 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
+	Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
+	Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2445,5 +2451,5 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
 	element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum);
-	Input* gllevelset_input = NULL;
+	Input2* gllevelset_input = NULL;
 
 	/*build friction object, used later on: */
@@ -2453,5 +2459,5 @@
 	if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base);
 	if(friction_style==SubelementFriction2Enum){
-		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		gllevelset_input=element->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
 		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
 		gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
@@ -2513,5 +2519,5 @@
 	/*Intermediaries*/
 	int         dim,bsize;
-	IssmDouble  viscosity,thickness,Jdet;
+	IssmDouble  viscosity,Jdet;
 	IssmDouble *xyz_list = NULL;
 
@@ -2529,8 +2535,8 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input    = element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input    = NULL;
+	Input2* vx_input    = element->GetInput2(VxEnum);       _assert_(vx_input);
+	Input2* vy_input    = NULL;
 	if(dim==3){
-		vy_input=element->GetInput(VyEnum);          _assert_(vy_input);
+		vy_input=element->GetInput2(VyEnum);          _assert_(vy_input);
 	}
 
@@ -2677,5 +2683,5 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input2*     surface_input = element->GetInput2(SurfaceEnum);   _assert_(surface_input);
 	IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum);
 
@@ -2733,6 +2739,6 @@
 
 	/*Retrieve all inputs and parameters*/
-	Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
+	Input2* sealevel_input       = element->GetInput2(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -2938,5 +2944,5 @@
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(PressureEnum,pressure,P1Enum);
 	xDelete<IssmDouble>(pressure);
 	xDelete<IssmDouble>(surface);
@@ -2983,7 +2989,7 @@
 
 	/*Add vx and vy as inputs to the element: */
-	element->AddInput(VxEnum,vx,element->GetElementType());
-	if(dim==3)element->AddInput(VyEnum,vy,element->GetElementType());
-	element->AddInput(VelEnum,vel,element->GetElementType());
+	element->AddInput2(VxEnum,vx,element->GetElementType());
+	if(dim==3)element->AddInput2(VyEnum,vy,element->GetElementType());
+	element->AddInput2(VelEnum,vel,element->GetElementType());
 
 	/*Free ressources:*/
@@ -3051,7 +3057,7 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
+	Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
+	Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
+	Input2* vz_input = element->GetInput2(VzEnum); _assert_(vz_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3141,5 +3147,5 @@
 	/*If on not water or not FS, skip stiffness: */
 	int approximation,shelf_dampening;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
 	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
@@ -3174,5 +3180,5 @@
 	IssmDouble  rho_water     = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble  gravity       = element->FindParam(ConstantsGEnum);
-	Input*      base_input = element->GetInput(BaseEnum); _assert_(base_input);
+	Input2*      base_input = element->GetInput2(BaseEnum); _assert_(base_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3232,8 +3238,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input = NULL;
-	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+	Input2* vx_input=element->GetInput2(VxEnum);     _assert_(vx_input);
+	Input2* vy_input=element->GetInput2(VyEnum);     _assert_(vy_input);
+	Input2* vz_input = NULL;
+	if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3382,8 +3388,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	//element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input = NULL;
-	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+	Input2* vx_input=element->GetInput2(VxEnum);     _assert_(vx_input);
+	Input2* vy_input=element->GetInput2(VyEnum);     _assert_(vy_input);
+	Input2* vz_input = NULL;
+	if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 
@@ -3456,8 +3462,8 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input = element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input = element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input = NULL;
-	if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
+	Input2* vx_input = element->GetInput2(VxEnum);     _assert_(vx_input);
+	Input2* vy_input = element->GetInput2(VyEnum);     _assert_(vy_input);
+	Input2* vz_input = NULL;
+	if(dim==3){vz_input = element->GetInput2(VzEnum); _assert_(vz_input);}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3581,8 +3587,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input;
-	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+	Input2* vx_input=element->GetInput2(VxEnum);     _assert_(vx_input);
+	Input2* vy_input=element->GetInput2(VyEnum);     _assert_(vy_input);
+	Input2* vz_input;
+	if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3623,5 +3629,5 @@
 	/*If on water or not FS, skip stiffness: */
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
 
@@ -3650,8 +3656,8 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
-	Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
-	Input* vz_input         = NULL;
-	if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
+	Input2* vx_input         = element->GetInput2(VxEnum);      _assert_(vx_input);
+	Input2* vy_input         = element->GetInput2(VyEnum);      _assert_(vy_input);
+	Input2* vz_input         = NULL;
+	if(dim==3){    vz_input = element->GetInput2(VzEnum);      _assert_(vz_input);}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3728,5 +3734,5 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	Input*  alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
+	Input2*  alpha2_input=element->GetInput2(FrictionCoefficientEnum); _assert_(alpha2_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3782,7 +3788,7 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	Input*  sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
-	Input*  sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
-	Input*  bedslope_input=element->GetInput(BedSlopeXEnum);     _assert_(bedslope_input);
+	Input2*  sigmann_input=element->GetInput2(VzEnum); _assert_(sigmann_input);
+	Input2*  sigmant_input=element->GetInput2(TemperatureEnum); _assert_(sigmant_input);
+	Input2*  bedslope_input=element->GetInput2(BedSlopeXEnum);     _assert_(bedslope_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3899,5 +3905,5 @@
 	/*If on water or not FS, skip stiffness: */
 	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
 
@@ -3927,5 +3933,5 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
 	element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum);
-	Input* gllevelset_input = NULL;
+	Input2* gllevelset_input = NULL;
 
 	/*build friction object, used later on: */
@@ -3936,5 +3942,5 @@
 	if(friction_style==SubelementFriction2Enum){
 		if(domaintype==Domain2DverticalEnum) _error_("Subelement Friction 2 not implemented yet for Flowline");
-		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		gllevelset_input=element->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
 		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
 		//gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
@@ -4068,9 +4074,9 @@
 	IssmDouble  rho_ice =element->FindParam(MaterialsRhoIceEnum);
 	IssmDouble  gravity =element->FindParam(ConstantsGEnum);
-	Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
-	Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
-	Input*      loadingforcez_input=NULL;
+	Input2*      loadingforcex_input=element->GetInput2(LoadingforceXEnum);  _assert_(loadingforcex_input);
+	Input2*      loadingforcey_input=element->GetInput2(LoadingforceYEnum);  _assert_(loadingforcey_input);
+	Input2*      loadingforcez_input=NULL;
 	if(dim==3){
-		loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
+		loadingforcez_input=element->GetInput2(LoadingforceZEnum);  _assert_(loadingforcez_input);
 	}
 
@@ -4146,6 +4152,6 @@
 	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
 	element->NormalSection(&normal[0],xyz_list_front);
-	Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input2* surface_input  = element->GetInput2(SurfaceEnum); _assert_(surface_input);
+	Input2* sealevel_input       = element->GetInput2(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble  rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble  gravity   = element->FindParam(ConstantsGEnum);
@@ -4215,5 +4221,5 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	Input*      base_input=element->GetInput(BaseEnum); _assert_(base_input);
+	Input2*      base_input=element->GetInput2(BaseEnum); _assert_(base_input);
 	IssmDouble  rho_water=element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble  gravity  =element->FindParam(ConstantsGEnum);
@@ -4239,5 +4245,5 @@
 	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
 	if(shelf_dampening) {
-		Input*      mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
+		Input2*      mb_input=element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
 		IssmDouble  dt,mb;
 		element->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -4290,6 +4296,6 @@
 
 	/*Get pressure and sigmann*/
-	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
+	Input2* pressure_input=element->GetInput2(PressureEnum); _assert_(pressure_input);
+	Input2* sigmann_input =element->GetInput2(SigmaNNEnum);  _assert_(sigmann_input);
 
 	gauss=element->NewGauss(5);
@@ -4383,19 +4389,19 @@
 
 	/*Get d and tau*/
-	Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
-	Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
-	Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
-	Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
-	Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
-	Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
-	Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
-	Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
+	Input2* epsxx_input=element->GetInput2(StrainRatexxEnum); _assert_(epsxx_input);
+	Input2* epsyy_input=element->GetInput2(StrainRateyyEnum); _assert_(epsyy_input);
+	Input2* epsxy_input=element->GetInput2(StrainRatexyEnum); _assert_(epsxy_input);
+	Input2* epszz_input=NULL; Input2* epsxz_input=NULL; Input2* epsyz_input=NULL;
+	Input2* sigmapxx_input=element->GetInput2(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
+	Input2* sigmapyy_input=element->GetInput2(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
+	Input2* sigmapxy_input=element->GetInput2(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
+	Input2* sigmapzz_input=NULL; Input2* sigmapxz_input=NULL; Input2* sigmapyz_input=NULL;
 	if(dim==3){
-		epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
-		epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
-		epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
-		sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
-		sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
-		sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
+		epszz_input=element->GetInput2(StrainRatezzEnum); _assert_(epszz_input);
+		epsxz_input=element->GetInput2(StrainRatexzEnum); _assert_(epsxz_input);
+		epsyz_input=element->GetInput2(StrainRateyzEnum); _assert_(epsyz_input);
+		sigmapzz_input=element->GetInput2(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
+		sigmapxz_input=element->GetInput2(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
+		sigmapyz_input=element->GetInput2(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
 	}
 
@@ -5059,5 +5065,5 @@
 	int*         vdoflist=NULL;
 	int*         pdoflist=NULL;
-	Input*       vz_input=NULL;
+	Input2*       vz_input=NULL;
 	int          dim;
 	IssmDouble   vx,vy,vz,p;
@@ -5081,8 +5087,8 @@
 	element->GetDofListVelocity(&vdoflist,GsetEnum);
 	element->GetDofListPressure(&pdoflist,GsetEnum);
-	Input*     vx_input=element->GetInput(VxEnum);       _assert_(vx_input);
-	Input*     vy_input=element->GetInput(VyEnum);       _assert_(vy_input);
-	if(dim==3){vz_input=element->GetInput(VzEnum);       _assert_(vz_input);}
-	Input*     p_input =element->GetInput(PressureEnum); _assert_(p_input);
+	Input2*     vx_input=element->GetInput2(VxEnum);       _assert_(vx_input);
+	Input2*     vy_input=element->GetInput2(VyEnum);       _assert_(vy_input);
+	if(dim==3){vz_input=element->GetInput2(VzEnum);       _assert_(vz_input);}
+	Input2*     p_input =element->GetInput2(PressureEnum); _assert_(p_input);
 
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
@@ -5142,8 +5148,8 @@
 		/*Get inputs and parameters*/
 		element->GetVerticesCoordinates(&xyz_list);
-		Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
-		Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
-		Input* vz_input;
-		if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+		Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
+		Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
+		Input2* vz_input;
+		if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 		/*Allocate new inputs*/
@@ -5184,11 +5190,11 @@
 
 		/*Add inputs*/
-		element->AddInput(StrainRatexxEnum,epsxx,P1DGEnum); element->AddInput(DeviatoricStressxxEnum,sigmapxx,P1DGEnum);
-		element->AddInput(StrainRateyyEnum,epsyy,P1DGEnum); element->AddInput(DeviatoricStressyyEnum,sigmapyy,P1DGEnum);
-		element->AddInput(StrainRatexyEnum,epsxy,P1DGEnum); element->AddInput(DeviatoricStressxyEnum,sigmapxy,P1DGEnum);
+		element->AddInput2(StrainRatexxEnum,epsxx,P1DGEnum); element->AddInput2(DeviatoricStressxxEnum,sigmapxx,P1DGEnum);
+		element->AddInput2(StrainRateyyEnum,epsyy,P1DGEnum); element->AddInput2(DeviatoricStressyyEnum,sigmapyy,P1DGEnum);
+		element->AddInput2(StrainRatexyEnum,epsxy,P1DGEnum); element->AddInput2(DeviatoricStressxyEnum,sigmapxy,P1DGEnum);
 		if(dim==3){
-			element->AddInput(StrainRatezzEnum,epszz,P1DGEnum); element->AddInput(DeviatoricStresszzEnum,sigmapzz,P1DGEnum);
-			element->AddInput(StrainRatexzEnum,epsxz,P1DGEnum); element->AddInput(DeviatoricStressxzEnum,sigmapxz,P1DGEnum);
-			element->AddInput(StrainRateyzEnum,epsyz,P1DGEnum); element->AddInput(DeviatoricStressyzEnum,sigmapyz,P1DGEnum);
+			element->AddInput2(StrainRatezzEnum,epszz,P1DGEnum); element->AddInput2(DeviatoricStresszzEnum,sigmapzz,P1DGEnum);
+			element->AddInput2(StrainRatexzEnum,epsxz,P1DGEnum); element->AddInput2(DeviatoricStressxzEnum,sigmapxz,P1DGEnum);
+			element->AddInput2(StrainRateyzEnum,epsyz,P1DGEnum); element->AddInput2(DeviatoricStressyzEnum,sigmapyz,P1DGEnum);
 		}
 
@@ -5273,12 +5279,12 @@
 
 	/*Add vx and vy as inputs to the tria element: */
-	int fe_v = element->VelocityInterpolation();
-	//if(fe_v == P1bubblecondensedEnum) fe_v = P1Enum;
-	//if(fe_v == P1bubbleEnum) fe_v = P1Enum;
-	element->AddInput(VxEnum, vx, fe_v);
-	element->AddInput(VyEnum, vy, fe_v);
-	element->AddInput(VelEnum,vel,fe_v);
-	if(pnumdof>0) element->AddInput(PressureEnum,pressure,element->PressureInterpolation());
-	if(dim==3) element->AddInput(VzEnum,vz, fe_v);
+	int v_interp =  element->VelocityInterpolation();
+	if(v_interp==P1bubbleEnum) v_interp=P1Enum;
+	if(v_interp == P1bubblecondensedEnum) v_interp = P1Enum;
+	element->AddInput2(VxEnum, vx, v_interp);
+	element->AddInput2(VyEnum, vy, v_interp);
+	element->AddInput2(VelEnum,vel,v_interp);
+	if(pnumdof>0) element->AddInput2(PressureEnum,pressure,element->PressureInterpolation());
+	if(dim==3) element->AddInput2(VzEnum,vz,v_interp);
 
 	/*Free ressources:*/
@@ -5314,10 +5320,10 @@
 		/*Get inputs and parameters*/
 		element->GetVerticesCoordinates(&xyz_list);
-		Input*  B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
-		Input*  n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
-		Input* vx_input=element->GetInput(VxEnum);                 _assert_(vx_input);
-		Input* vy_input=element->GetInput(VyEnum);                 _assert_(vy_input);
-		Input* vz_input;
-		if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+		Input2*  B_input=element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
+		Input2*  n_input=element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
+		Input2* vx_input=element->GetInput2(VxEnum);                 _assert_(vx_input);
+		Input2* vy_input=element->GetInput2(VyEnum);                 _assert_(vy_input);
+		Input2* vz_input;
+		if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 		/*Fetch number of nodes and dof for this finite element*/
@@ -5338,23 +5344,23 @@
 
 		/*Get previous d*/
-		Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
-		Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
-		Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
-		Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+		Input2* epsxx_input=element->GetInput2(StrainRatexxEnum); _assert_(epsxx_input);
+		Input2* epsyy_input=element->GetInput2(StrainRateyyEnum); _assert_(epsyy_input);
+		Input2* epsxy_input=element->GetInput2(StrainRatexyEnum); _assert_(epsxy_input);
+		Input2* epszz_input=NULL; Input2* epsxz_input=NULL; Input2* epsyz_input=NULL;
 		if(dim==3){
-			epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
-			epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
-			epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+			epszz_input=element->GetInput2(StrainRatezzEnum); _assert_(epszz_input);
+			epsxz_input=element->GetInput2(StrainRatexzEnum); _assert_(epsxz_input);
+			epsyz_input=element->GetInput2(StrainRateyzEnum); _assert_(epsyz_input);
 		}
 
 		/*Get tau*/
-		Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
-		Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
-		Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
-		Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
+		Input2* sigmapxx_input=element->GetInput2(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
+		Input2* sigmapyy_input=element->GetInput2(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
+		Input2* sigmapxy_input=element->GetInput2(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
+		Input2* sigmapzz_input=NULL; Input2* sigmapxz_input=NULL; Input2* sigmapyz_input=NULL;
 		if(dim==3){
-			sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
-			sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
-			sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
+			sigmapzz_input=element->GetInput2(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
+			sigmapxz_input=element->GetInput2(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
+			sigmapyz_input=element->GetInput2(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
 		}
 
@@ -5470,7 +5476,7 @@
 			for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_yy[i]));
 			for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
-			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
-			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
-			element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
+			element->AddInput2(StrainRatexxEnum,d_xx,P1DGEnum);
+			element->AddInput2(StrainRateyyEnum,d_yy,P1DGEnum);
+			element->AddInput2(StrainRatexyEnum,d_xy,P1DGEnum);
 		}
 		else{
@@ -5485,10 +5491,10 @@
 			Matrix4x4Solve(&d_xz[0],Ke,pe_xz);
 			Matrix4x4Solve(&d_yz[0],Ke,pe_yz);
-			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
-			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
-			element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
-			element->AddInput(StrainRatezzEnum,d_zz,P1DGEnum);
-			element->AddInput(StrainRatexzEnum,d_xz,P1DGEnum);
-			element->AddInput(StrainRateyzEnum,d_yz,P1DGEnum);
+			element->AddInput2(StrainRatexxEnum,d_xx,P1DGEnum);
+			element->AddInput2(StrainRateyyEnum,d_yy,P1DGEnum);
+			element->AddInput2(StrainRatexyEnum,d_xy,P1DGEnum);
+			element->AddInput2(StrainRatezzEnum,d_zz,P1DGEnum);
+			element->AddInput2(StrainRatexzEnum,d_xz,P1DGEnum);
+			element->AddInput2(StrainRateyzEnum,d_yz,P1DGEnum);
 		}
 
@@ -5527,29 +5533,29 @@
 		/*Get inputs and parameters*/
 		element->GetVerticesCoordinates(&xyz_list);
-		Input* vx_input=element->GetInput(VxEnum);                 _assert_(vx_input);
-		Input* vy_input=element->GetInput(VyEnum);                 _assert_(vy_input);
-		Input* vz_input=NULL;
-		if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+		Input2* vx_input=element->GetInput2(VxEnum);                 _assert_(vx_input);
+		Input2* vy_input=element->GetInput2(VyEnum);                 _assert_(vy_input);
+		Input2* vz_input=NULL;
+		if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
 
 		/*Get previous tau*/
-		Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
-		Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
-		Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
-		Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
+		Input2* sigmapxx_input=element->GetInput2(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
+		Input2* sigmapyy_input=element->GetInput2(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
+		Input2* sigmapxy_input=element->GetInput2(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
+		Input2* sigmapzz_input=NULL; Input2* sigmapxz_input=NULL; Input2* sigmapyz_input=NULL;
 		if(dim==3){
-			sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
-			sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
-			sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
+			sigmapzz_input=element->GetInput2(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
+			sigmapxz_input=element->GetInput2(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
+			sigmapyz_input=element->GetInput2(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
 		}
 
 		/*Get NEW d*/
-		Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
-		Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
-		Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
-		Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+		Input2* epsxx_input=element->GetInput2(StrainRatexxEnum); _assert_(epsxx_input);
+		Input2* epsyy_input=element->GetInput2(StrainRateyyEnum); _assert_(epsyy_input);
+		Input2* epsxy_input=element->GetInput2(StrainRatexyEnum); _assert_(epsxy_input);
+		Input2* epszz_input=NULL; Input2* epsxz_input=NULL; Input2* epsyz_input=NULL;
 		if(dim==3){
-			epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
-			epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
-			epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+			epszz_input=element->GetInput2(StrainRatezzEnum); _assert_(epszz_input);
+			epsxz_input=element->GetInput2(StrainRatexzEnum); _assert_(epsxz_input);
+			epsyz_input=element->GetInput2(StrainRateyzEnum); _assert_(epsyz_input);
 		}
 
@@ -5620,11 +5626,11 @@
 
 		/*Add inputs*/
-		element->AddInput(DeviatoricStressxxEnum,tau_xx,P1DGEnum);
-		element->AddInput(DeviatoricStressyyEnum,tau_yy,P1DGEnum);
-		element->AddInput(DeviatoricStressxyEnum,tau_xy,P1DGEnum);
+		element->AddInput2(DeviatoricStressxxEnum,tau_xx,P1DGEnum);
+		element->AddInput2(DeviatoricStressyyEnum,tau_yy,P1DGEnum);
+		element->AddInput2(DeviatoricStressxyEnum,tau_xy,P1DGEnum);
 		if(dim==3){
-			element->AddInput(DeviatoricStresszzEnum,tau_zz,P1DGEnum);
-			element->AddInput(DeviatoricStressxzEnum,tau_xz,P1DGEnum);
-			element->AddInput(DeviatoricStressyzEnum,tau_yz,P1DGEnum);
+			element->AddInput2(DeviatoricStresszzEnum,tau_zz,P1DGEnum);
+			element->AddInput2(DeviatoricStressxzEnum,tau_xz,P1DGEnum);
+			element->AddInput2(DeviatoricStressyzEnum,tau_yz,P1DGEnum);
 		}
 
@@ -5752,5 +5758,5 @@
 
 	/*If on water or not FS, skip stiffness: */
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(element->IsFloating() || !element->IsOnBase()) return NULL;
 
@@ -5784,7 +5790,7 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_tria);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
+	Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
+	Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
+	Input2* vz_input=element->GetInput2(VzEnum); _assert_(vz_input);
 
 	/*build friction object, used later on: */
@@ -5905,7 +5911,7 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
+	Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
+	Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
+	Input2* vz_input=element->GetInput2(VzEnum); _assert_(vz_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -6092,6 +6098,6 @@
 	/* Get node coordinates and dof list: */
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
+	Input2* vx_input   =element->GetInput2(VxEnum);       _assert_(vx_input);
+	Input2* vy_input   =element->GetInput2(VyEnum);       _assert_(vy_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -6245,11 +6251,11 @@
 	/*Initialize Element matrix*/
 	ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
-	Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
+	Input2* vx_input   =element->GetInput2(VxEnum);       _assert_(vx_input);
+	Input2* vy_input   =element->GetInput2(VyEnum);       _assert_(vy_input);
+	Input2* vz_input   =element->GetInput2(VzEnum);       _assert_(vz_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -6326,5 +6332,5 @@
 	/*Initialize Element vector and return if necessary*/
 	if(!element->IsOnBase() || element->IsFloating()) return NULL;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=HOFSApproximationEnum) return NULL;
 
@@ -6351,8 +6357,8 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_tria);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=  element->GetInput(VxEnum);   _assert_(vx_input);
-	Input* vy_input=  element->GetInput(VyEnum);   _assert_(vy_input);
-	Input* vz_input=  element->GetInput(VzEnum);   _assert_(vz_input);
-	Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input);
+	Input2* vx_input=  element->GetInput2(VxEnum);   _assert_(vx_input);
+	Input2* vy_input=  element->GetInput2(VyEnum);   _assert_(vy_input);
+	Input2* vz_input=  element->GetInput2(VzEnum);   _assert_(vz_input);
+	Input2* vzHO_input=element->GetInput2(VzHOEnum); _assert_(vzHO_input);
 
 	/*build friction object, used later on: */
@@ -6409,5 +6415,5 @@
 
 	/*Initialize Element vector and return if necessary*/
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=HOFSApproximationEnum) return NULL;
 	int   vnumnodes = element->NumberofNodesVelocity();
@@ -6430,8 +6436,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input   =element->GetInput(VxEnum);   _assert_(vx_input);
-	Input* vy_input   =element->GetInput(VyEnum);   _assert_(vy_input);
-	Input* vz_input   =element->GetInput(VzEnum);   _assert_(vz_input);
-	Input* vzHO_input=element->GetInput(VzHOEnum);  _assert_(vzHO_input);
+	Input2* vx_input   =element->GetInput2(VxEnum);   _assert_(vx_input);
+	Input2* vy_input   =element->GetInput2(VyEnum);   _assert_(vy_input);
+	Input2* vz_input   =element->GetInput2(VzEnum);   _assert_(vz_input);
+	Input2* vzHO_input=element->GetInput2(VzHOEnum);  _assert_(vzHO_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -6496,5 +6502,5 @@
 	/*Initialize Element vector and return if necessary*/
 	if(!element->IsOnBase() || element->IsFloating()) return NULL;
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=SSAFSApproximationEnum) return NULL;
 	int vnumnodes = element->NumberofNodesVelocity();
@@ -6518,8 +6524,8 @@
 	element->GetVerticesCoordinatesBase(&xyz_list_tria);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input=   element->GetInput(VxEnum);    _assert_(vx_input);
-	Input* vy_input=   element->GetInput(VyEnum);    _assert_(vy_input);
-	Input* vz_input=   element->GetInput(VzEnum);    _assert_(vz_input);
-	Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input);
+	Input2* vx_input=   element->GetInput2(VxEnum);    _assert_(vx_input);
+	Input2* vy_input=   element->GetInput2(VyEnum);    _assert_(vy_input);
+	Input2* vz_input=   element->GetInput2(VzEnum);    _assert_(vz_input);
+	Input2* vzSSA_input=element->GetInput2(VzSSAEnum); _assert_(vzSSA_input);
 
 	/*build friction object, used later on: */
@@ -6575,5 +6581,5 @@
 
 	/*Initialize Element vector and return if necessary*/
-	element->GetInputValue(&approximation,ApproximationEnum);
+	element->GetInput2Value(&approximation,ApproximationEnum);
 	if(approximation!=SSAFSApproximationEnum) return NULL;
 	int vnumnodes = element->NumberofNodesVelocity();
@@ -6596,8 +6602,8 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	Input* vx_input   =element->GetInput(VxEnum);      _assert_(vx_input);
-	Input* vy_input   =element->GetInput(VyEnum);      _assert_(vy_input);
-	Input* vz_input   =element->GetInput(VzEnum);      _assert_(vz_input);
-	Input* vzSSA_input=element->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
+	Input2* vx_input   =element->GetInput2(VxEnum);      _assert_(vx_input);
+	Input2* vy_input   =element->GetInput2(VyEnum);      _assert_(vy_input);
+	Input2* vz_input   =element->GetInput2(VzEnum);      _assert_(vz_input);
+	Input2* vzSSA_input=element->GetInput2(VzSSAEnum);   _assert_(vzSSA_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -7192,10 +7198,10 @@
 
 	/*Add vx and vy as inputs to element: */
-	element->AddInput(VxEnum,vx,P1Enum);
-	element->AddInput(VyEnum,vy,P1Enum);
-	element->AddInput(VzEnum,vz,P1Enum);
-	element->AddInput(VzFSEnum,vzFS,P1Enum);
-	element->AddInput(VelEnum,vel,P1Enum);
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(VxEnum,vx,P1Enum);
+	element->AddInput2(VyEnum,vy,P1Enum);
+	element->AddInput2(VzEnum,vz,P1Enum);
+	element->AddInput2(VzFSEnum,vzFS,P1Enum);
+	element->AddInput2(VelEnum,vel,P1Enum);
+	//element->AddInput2(PressureEnum,pressure,P1Enum);
 
 	/*Free ressources:*/
@@ -7295,10 +7301,10 @@
 
 	/*Add vx and vy as inputs to element: */
-	element->AddInput(VxEnum,vx,P1Enum);
-	element->AddInput(VyEnum,vy,P1Enum);
-	element->AddInput(VzEnum,vz,P1Enum);
-	element->AddInput(VzFSEnum,vzFS,P1Enum);
-	element->AddInput(VelEnum,vel,P1Enum);
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(VxEnum,vx,P1Enum);
+	element->AddInput2(VyEnum,vy,P1Enum);
+	element->AddInput2(VzEnum,vz,P1Enum);
+	element->AddInput2(VzFSEnum,vzFS,P1Enum);
+	element->AddInput2(VelEnum,vel,P1Enum);
+	//element->AddInput2(PressureEnum,pressure,P1Enum);
 
 	/*Free ressources:*/
@@ -7388,8 +7394,8 @@
 
 	/*Add vx and vy as inputs to element: */
-	element->AddInput(VxEnum,vx,P1Enum);
-	element->AddInput(VyEnum,vy,P1Enum);
-	element->AddInput(VelEnum,vel,P1Enum);
-	element->AddInput(PressureEnum,pressure,P1Enum);
+	element->AddInput2(VxEnum,vx,P1Enum);
+	element->AddInput2(VyEnum,vy,P1Enum);
+	element->AddInput2(VelEnum,vel,P1Enum);
+	element->AddInput2(PressureEnum,pressure,P1Enum);
 
 	/*Free ressources:*/
