Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17275)
@@ -270,6 +270,6 @@
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 	this->CreateDamageFInput(element);
-	Input* damaged_input = element->GetMaterialInput(DamageDbarEnum); _assert_(damaged_input);
-	Input* damagef_input = element->GetMaterialInput(DamageFEnum);    _assert_(damagef_input);
+	Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
+	Input* damagef_input = element->GetInput(DamageFEnum);    _assert_(damagef_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -378,5 +378,5 @@
 
 	/*Get all inputs and parameters*/
-	element->AddMaterialInput(DamageDbarEnum,values,P1Enum);
+	element->AddInput(DamageDbarEnum,values,P1Enum);
 
 	/*Free ressources:*/
@@ -416,5 +416,5 @@
 	Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
 	Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
-	Input* damage_input    = element->GetMaterialInput(DamageDbarEnum); _assert_(damage_input);
+	Input* damage_input    = element->GetInput(DamageDbarEnum); _assert_(damage_input);
 
 	/*Damage evolution z mapping: */
@@ -442,5 +442,5 @@
 
 	/*Add input*/
-	element->AddMaterialInput(DamageFEnum,f,P1Enum);
+	element->AddInput(DamageFEnum,f,P1Enum);
 	
 	/*Clean up and return*/
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17275)
@@ -790,14 +790,14 @@
 			case PatersonEnum:
 				for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
+				element->AddInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 			case ArrheniusEnum:
 				element->GetVerticesCoordinates(&xyz_list);
 				for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
+				element->AddInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 			case LliboutryDuvalEnum:
 				for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],element->GetMaterialParameter(MaterialsRheologyNEnum),element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum)); 
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum); 
+				element->AddInput(MaterialsRheologyBEnum,&B[0],P1Enum); 
 				break; 
 			default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17275)
@@ -1044,5 +1044,5 @@
 void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
-	_printf0_("   Updating active and non-active nodes for StressbalanceAnalysis \n");
+	//_printf0_("   Updating active and non-active nodes for StressbalanceAnalysis \n");
 //	SetActiveNodesLSMx(femmodel->elements);
 	return;
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 17275)
@@ -659,10 +659,10 @@
 			case PatersonEnum:
 				for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
+				element->AddInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 			case ArrheniusEnum:{
 				element->GetVerticesCoordinates(&xyz_list);
 				for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
+				element->AddInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 				}
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17275)
@@ -155,6 +155,5 @@
 			return this->matpar->GetMaterialParameter(enum_in);
 	}
-}
-/*}}}*/
+}/*}}}*/
 void Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
 	/*Compute deformational heating from epsilon and viscosity */
@@ -227,4 +226,79 @@
 	return shelf;
 }/*}}}*/
+bool Element::IsInput(int name){/*{{{*/
+	if (
+				name==ThicknessEnum ||
+				name==SurfaceEnum ||
+				name==BedEnum ||
+				name==BathymetryEnum ||
+				name==SurfaceSlopeXEnum ||
+				name==SurfaceSlopeYEnum ||
+				name==SurfaceforcingsMassBalanceEnum ||
+				name==BasalforcingsMeltingRateEnum ||
+				name==BasalforcingsGeothermalfluxEnum ||
+				name==SurfaceAreaEnum||
+				name==DamageDEnum ||
+				name==DamageDbarEnum ||
+				name==PressureEnum ||
+				name==VxEnum ||
+				name==VyEnum ||
+				name==VzEnum ||
+				name==VxMeshEnum ||
+				name==VyMeshEnum ||
+				name==VzMeshEnum ||
+				name==InversionVxObsEnum ||
+				name==InversionVyObsEnum ||
+				name==InversionVzObsEnum ||
+				name==TemperatureEnum ||
+				name==EnthalpyEnum ||
+				name==EnthalpyPicardEnum ||
+				name==WaterfractionEnum||
+				name==WatercolumnEnum || 
+				name==FrictionCoefficientEnum ||
+				name==MaskGroundediceLevelsetEnum ||
+				name==MaskIceLevelsetEnum ||
+				name==IceMaskNodeActivationEnum ||
+				name==LevelsetfunctionSlopeXEnum ||
+				name==LevelsetfunctionSlopeYEnum ||
+				name==GradientEnum ||
+				name==OldGradientEnum  ||
+				name==ConvergedEnum || 
+				name==QmuVxEnum ||
+				name==QmuVyEnum ||
+				name==QmuVzEnum ||
+				name==QmuVxMeshEnum ||
+				name==QmuVyMeshEnum ||
+				name==QmuVzMeshEnum ||
+				name==QmuPressureEnum ||
+				name==QmuBedEnum ||
+				name==QmuThicknessEnum ||
+				name==QmuSurfaceEnum ||
+				name==QmuTemperatureEnum ||
+				name==QmuMeltingEnum ||
+				name==QmuMaskGroundediceLevelsetEnum ||
+				name==QmuMaskIceLevelsetEnum ||
+				name==QmuMaterialsRheologyBEnum ||
+				name==MaterialsRheologyBEnum ||
+				name==MaterialsRheologyBbarEnum ||
+				name==MaterialsRheologyNEnum ||
+				name==GiaWEnum || 
+				name==GiadWdtEnum ||
+				name==SedimentHeadEnum ||
+				name==EplHeadEnum ||
+				name==SedimentHeadOldEnum ||
+				name==EplHeadOldEnum ||
+				name==HydrologydcEplThicknessOldEnum ||
+				name==HydrologydcEplInitialThicknessEnum ||
+				name==HydrologydcEplThicknessEnum ||
+				name==HydrologydcMaskEplactiveNodeEnum ||
+				name==MeshVertexonbedEnum ||
+				name==WaterTransferEnum
+
+				) {
+					return true;
+				}
+	else return false;
+}
+/*}}}*/
 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17275)
@@ -65,4 +65,5 @@
 		void       GetVerticesSidList(int* sidlist);
 		void       GetVerticesConnectivityList(int* connectivitylist);
+		bool	     IsInput(int name);
 		bool       IsFloating(); 
 		void       ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
@@ -106,5 +107,4 @@
 		virtual void       AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual void       AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
-		virtual void       AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual IssmDouble CharacteristicLength(void)=0;
 		virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
@@ -164,5 +164,4 @@
 		virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype)=0;
 		virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
-		virtual Input* GetMaterialInput(int inputenum)=0;
 		virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 		virtual void   GetInputValue(bool* pvalue,int enum_type)=0;
@@ -177,9 +176,8 @@
 		virtual IssmDouble GetYcoord(Gauss* gauss)=0;
 		virtual IssmDouble GetZcoord(Gauss* gauss)=0;
-		virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 		virtual int    GetElementType(void)=0;
 
 		virtual IssmDouble SurfaceArea(void)=0;
-		virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;
+		virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
 		virtual void   InputChangeName(int enum_type,int enum_type_old)=0;
 		virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17275)
@@ -139,11 +139,4 @@
 		else _error_("not implemented yet");
 	}
-}
-/*}}}*/
-/*FUNCTION Penta::AddMaterialInput{{{*/
-void  Penta::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
-
-	_assert_(this->material);
-	this->material->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
@@ -465,6 +458,6 @@
 	this->inputs->AddInput(NewPrecipitationInput);
 
-	this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,ElementEnum);
-	this->InputExtrude(SurfaceforcingsPrecipitationEnum,ElementEnum);
+	this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum);
+	this->InputExtrude(SurfaceforcingsPrecipitationEnum);
 
 	/*clean-up*/
@@ -974,9 +967,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetMaterialInput(int inputenum) {{{*/
-Input* Penta::GetMaterialInput(int inputenum){
-	return this->material->inputs->GetInput(inputenum);
-}
-/*}}}*/
 /*FUNCTION Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype) {{{*/
 void Penta::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){
@@ -1138,17 +1126,4 @@
 
 }/*}}}*/
-/*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
-void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
-
-	Input* input=this->material->inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
-
-	GaussPenta* gauss=new GaussPenta();
-	gauss->GaussVertex(this->GetNodeIndex(node));
-
-	input->GetInputValue(pvalue,gauss);
-	delete gauss;
-}
-/*}}}*/
 /*FUNCTION Penta::NormalSection{{{*/
 void Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
@@ -1500,5 +1475,5 @@
 /*}}}*/
 /*FUNCTION Penta::InputDepthAverageAtBase{{{*/
-void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
+void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
 
 	int  step,i;
@@ -1528,10 +1503,5 @@
 
 		/*Step1: Get original input (to be depth avegaged): */
-		if (object_enum==MeshElementsEnum)
-		 original_input=(Input*)penta->inputs->GetInput(enum_type);
-		else if (object_enum==MaterialsEnum)
-		 original_input=(Input*)penta->material->inputs->GetInput(enum_type);
-		else
-		 _error_("object " << EnumToStringx(object_enum) << " not supported yet");
+		original_input=(Input*)penta->inputs->GetInput(enum_type);
 		if(!original_input) _error_("could not find input with enum " << EnumToStringx(enum_type));
 
@@ -1589,10 +1559,5 @@
 
 	/*Finally, add to inputs*/
-	if (object_enum==MeshElementsEnum)
-	 this->inputs->AddInput((Input*)depth_averaged_input);
-	else if (object_enum==MaterialsEnum)
-	 this->material->inputs->AddInput((Input*)depth_averaged_input);
-	else
-	 _error_("object " << EnumToStringx(object_enum) << " not supported yet");
+	this->inputs->AddInput((Input*)depth_averaged_input);
 }
 /*}}}*/
@@ -1606,5 +1571,5 @@
 /*}}}*/
 /*FUNCTION Penta::InputExtrude {{{*/
-void  Penta::InputExtrude(int enum_type,int object_type){
+void  Penta::InputExtrude(int enum_type){
 
 	int     i,num_inputs;
@@ -1617,19 +1582,9 @@
 
 	/*Step1: Get and Extrude original input: */
-	if(object_type==ElementEnum){
-		num_inputs=1;
-		base_inputs=xNew<Input*>(num_inputs);
-		base_inputs[0]=(Input*)this->inputs->GetInput(enum_type);
-	}
-	else if(object_type==MaterialsEnum){
-		num_inputs=1;
-		base_inputs=xNew<Input*>(num_inputs);
-		base_inputs[0]=(Input*)material->inputs->GetInput(enum_type);
-	}
-	else{
-		_error_("object of type " << EnumToStringx(object_type) << " not supported yet");
-	}
+	num_inputs=1;
+	base_inputs=xNew<Input*>(num_inputs);
+	base_inputs[0]=(Input*)this->inputs->GetInput(enum_type);
 	for(i=0;i<num_inputs;i++){
-		if(!base_inputs[i]) _error_("could not find input with enum " << EnumToStringx(enum_type) << " in object " << EnumToStringx(object_type));
+		if(!base_inputs[i]) _error_("could not find input with enum " << EnumToStringx(enum_type));
 		base_inputs[i]->Extrude();
 	}
@@ -1648,13 +1603,5 @@
 		for(i=0;i<num_inputs;i++){
 			copy=(Input*)base_inputs[i]->copy();
-			if (object_type==ElementEnum){
-				penta->inputs->AddInput((Input*)copy);
-			}
-			else if(object_type==MaterialsEnum){
-				penta->material->inputs->AddInput((Input*)copy);
-			}
-			else{
-				_error_("object of type " << EnumToStringx(object_type) << " not supported yet");
-			}
+			penta->inputs->AddInput((Input*)copy);
 		}
 
@@ -1774,7 +1721,20 @@
 					}
 					break;
-				/*Material will take care of it*/ 
-				case MaterialsRheologyBbarEnum: break;
-				case DamageDbarEnum:break;
+				case MaterialsRheologyBbarEnum:
+					if(iomodel->Data(MaterialsRheologyBEnum)){
+						for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1];
+						for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+						for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+						this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+					}
+					break;
+				case DamageDbarEnum:
+					if(iomodel->Data(DamageDEnum)){
+						for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[penta_vertex_ids[j]-1];
+						for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+						for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
+						this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+					}
+					break;
 				default:
 					_error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
@@ -1944,72 +1904,4 @@
 	inputs->GetInputValue(&onsurface,MeshElementonsurfaceEnum);
 	return onsurface;
-}
-/*}}}*/
-/*FUNCTION Penta::IsInput{{{*/
-bool Penta::IsInput(int name){
-	if (
-				name==ThicknessEnum ||
-				name==SurfaceEnum ||
-				name==BedEnum ||
-				name==BathymetryEnum ||
-				name==SurfaceSlopeXEnum ||
-				name==SurfaceSlopeYEnum ||
-				name==SurfaceforcingsMassBalanceEnum ||
-				name==BasalforcingsMeltingRateEnum ||
-				name==BasalforcingsGeothermalfluxEnum ||
-				name==SurfaceAreaEnum||
-				name==PressureEnum ||
-				name==VxEnum ||
-				name==VyEnum ||
-				name==VzEnum ||
-				name==VxMeshEnum ||
-				name==VyMeshEnum ||
-				name==VzMeshEnum ||
-				name==InversionVxObsEnum ||
-				name==InversionVyObsEnum ||
-				name==InversionVzObsEnum ||
-				name==TemperatureEnum ||
-				name==EnthalpyEnum ||
-				name==EnthalpyPicardEnum ||
-				name==WaterfractionEnum||
-				name==FrictionCoefficientEnum ||
-				name==MaskGroundediceLevelsetEnum ||
-				name==MaskIceLevelsetEnum ||
-				name==IceMaskNodeActivationEnum ||
-				name==LevelsetfunctionSlopeXEnum ||
-				name==LevelsetfunctionSlopeYEnum ||
-				name==GradientEnum ||
-				name==OldGradientEnum  ||
-				name==ConvergedEnum || 
-				name==QmuVxEnum ||
-				name==QmuVyEnum ||
-				name==QmuVzEnum ||
-				name==QmuVxMeshEnum ||
-				name==QmuVyMeshEnum ||
-				name==QmuVzMeshEnum ||
-				name==QmuPressureEnum ||
-				name==QmuBedEnum ||
-				name==QmuThicknessEnum ||
-				name==QmuSurfaceEnum ||
-				name==QmuTemperatureEnum ||
-				name==QmuMeltingEnum ||
-				name==QmuMaskGroundediceLevelsetEnum ||
-				name==QmuMaskIceLevelsetEnum ||
-				name==GiaWEnum || 
-				name==GiadWdtEnum ||
-				name==SedimentHeadEnum ||
-				name==EplHeadEnum ||
-				name==SedimentHeadOldEnum ||
-				name==EplHeadOldEnum ||
-				name==HydrologydcEplThicknessOldEnum ||
-				name==HydrologydcEplInitialThicknessEnum ||
-				name==HydrologydcEplThicknessEnum ||
-				name==HydrologydcMaskEplactiveNodeEnum ||
-				name==WaterTransferEnum
-
-				) {
-		return true;
-	}
-	else return false;
 }
 /*}}}*/
@@ -2354,5 +2246,5 @@
    this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
-   this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
+   this->InputExtrude(SurfaceforcingsMassBalanceEnum);
 
 	/*clean-up*/
@@ -2579,7 +2471,5 @@
 
 	/*Spawn material*/
-	tria->material=(Material*)this->material->copy();
-	delete tria->material->inputs;
-	tria->material->inputs=(Inputs*)this->material->inputs->SpawnTriaInputs(location);
+	tria->material=(Material*)this->material->copy(tria);
 
 	/*recover nodes, material and matpar: */
@@ -2597,11 +2487,11 @@
 	_assert_(this->IsOnBed());
 
-	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
-	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 	if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
 	if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 	Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
-	this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
-	this->material->inputs->DeleteInput(DamageDbarEnum);
+	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+	this->inputs->DeleteInput(DamageDbarEnum);
 	this->inputs->DeleteInput(VxAverageEnum);
 	this->inputs->DeleteInput(VyAverageEnum);
@@ -3491,9 +3381,9 @@
 	if(enum_type==MaterialsRheologyBbarEnum){
 		if(!IsOnBed()) return;
-		input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
 	}
 	else if(enum_type==DamageDbarEnum){
 		if(!IsOnBed()) return;
-		input=(Input*)material->inputs->GetInput(DamageDEnum);
+		input=(Input*)inputs->GetInput(DamageDEnum);
 	}
 	else{
@@ -3513,8 +3403,8 @@
 
 	if(enum_type==MaterialsRheologyBbarEnum){
-		input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
 	}
 	else if(enum_type==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(DamageDEnum);
+		input=(Input*)inputs->GetInput(DamageDEnum);
 	}
 	else{
@@ -3535,8 +3425,8 @@
 
 	if(enum_type==MaterialsRheologyBbarEnum){
-		input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
 	}
 	else if(enum_type==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(DamageDEnum);
+		input=(Input*)inputs->GetInput(DamageDEnum);
 	}
 	else{
@@ -3558,8 +3448,8 @@
 
 	if(control_enum==MaterialsRheologyBbarEnum){
-		input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
 	}
 	else if(control_enum==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(DamageDEnum);
+		input=(Input*)inputs->GetInput(DamageDEnum);
 	}
 	else{
@@ -3883,6 +3773,6 @@
 
 	/*Depth Average B*/
-	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
-	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 
 	/*Collapse element to the base*/
@@ -3892,6 +3782,6 @@
 
 	/*delete Average B*/
-	this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
-	this->material->inputs->DeleteInput(DamageDbarEnum);
+	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+	this->inputs->DeleteInput(DamageDbarEnum);
 
 } /*}}}*/
@@ -3903,6 +3793,6 @@
 
 	/*Depth Average B and D*/
-	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
-	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 
 	/*Collapse element to the base*/
@@ -3912,6 +3802,6 @@
 
 	/*delete Average B*/
-	this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
-	this->material->inputs->DeleteInput(DamageDbarEnum);
+	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+	this->inputs->DeleteInput(DamageDbarEnum);
 } /*}}}*/
 /*FUNCTION Penta::GradjBbarFS {{{*/
@@ -3922,6 +3812,6 @@
 
 	/*Depth Average B and D*/
-	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
-	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 
 	/*Collapse element to the base*/
@@ -3931,6 +3821,6 @@
 
 	/*delete Average B*/
-	this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
-	this->material->inputs->DeleteInput(DamageDbarEnum);
+	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+	this->inputs->DeleteInput(DamageDbarEnum);
 } /*}}}*/
 /*FUNCTION Penta::InputControlUpdate{{{*/
@@ -3950,15 +3840,14 @@
 		if(control_type[i]==MaterialsRheologyBbarEnum){
 			if (!IsOnBed()) goto cleanup_and_return;
-			input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
+			input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
 		}
 		else if(control_type[i]==DamageDbarEnum){
 			if (!IsOnBed()) goto cleanup_and_return;
-			input=(Input*)material->inputs->GetInput(DamageDEnum); _assert_(input);
+			input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
 		}
 		else{
 			input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
 		}
-
-		if (input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
+		if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
 
 		((ControlInput*)input)->UpdateValue(scalar);
@@ -3967,8 +3856,8 @@
 
 		if(control_type[i]==MaterialsRheologyBbarEnum){
-			this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum);
+			this->InputExtrude(MaterialsRheologyBEnum);
 		}
 		else if(control_type[i]==DamageDbarEnum){
-			this->InputExtrude(DamageDEnum,MaterialsEnum);
+			this->InputExtrude(DamageDEnum);
 		}
 	}
@@ -4264,11 +4153,5 @@
 	}
 	new_input = new PentaInput(control_enum,values,P1Enum);
-
-	if(control_enum==MaterialsRheologyBbarEnum){
-		input=(Input*)material->inputs->GetInput(control_enum); _assert_(input);
-	}
-	else{
-		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
-	}
+	input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
 
 	if (input->ObjectEnum()!=ControlInputEnum){
@@ -4369,6 +4252,6 @@
 					this->inputs->AddInput(new PentaInput(BedEnum,bed,P1Enum));
 					this->inputs->AddInput(new PentaInput(SurfaceEnum,surface,P1Enum));
-
 					break;
+
 				default:
 					this->inputs->AddInput(new PentaInput(name,values,P1Enum));
@@ -4656,5 +4539,5 @@
 	}
 	this->inputs->AddInput(new PentaInput(MaskGroundediceLevelsetEnum,&phi[0],P1Enum));
-	this->InputExtrude(MaskGroundediceLevelsetEnum,ElementEnum);
+	this->InputExtrude(MaskGroundediceLevelsetEnum);
 
 	/*SubelementMigrationEnum: if one grounded, all grounded*/
@@ -4688,6 +4571,6 @@
 
 	/*Extrude inputs*/
-	this->InputExtrude(SurfaceEnum,ElementEnum);
-	this->InputExtrude(BedEnum,ElementEnum);
+	this->InputExtrude(SurfaceEnum);
+	this->InputExtrude(BedEnum);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17275)
@@ -105,5 +105,5 @@
 		int    Sid();
 		void   InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
-		void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
+		void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
 		void   InputDuplicate(int original_enum,int new_enum);
 		void   InputScale(int enum_type,IssmDouble scale_factor);
@@ -195,5 +195,4 @@
 		void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
-		void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		void	         NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
 		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
@@ -206,5 +205,4 @@
 		void	         GetVertexPidList(int* doflist);
 		int            GetElementType(void);
-		Input*         GetMaterialInput(int inputenum);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
@@ -216,11 +214,9 @@
 		void           GetInputValue(IssmDouble* pvalue,int enum_type);
 		void           GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
-		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		Node*          GetNode(int node_number);
 		void           InputChangeName(int input_enum, int enum_type_old);
-		void	         InputExtrude(int enum_type,int object_type);
+		void	         InputExtrude(int enum_type);
 		void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
 		void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
-		bool	         IsInput(int name);
 		bool	         IsOnSurface(void);
 		bool	         IsOnBed(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17275)
@@ -61,5 +61,4 @@
 		void        AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
 		void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
-		void        AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
 		IssmDouble  CharacteristicLength(void);
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
@@ -130,5 +129,4 @@
 		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
 		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
-		Input*      GetMaterialInput(int inputenum){_error_("not implemented yet");};
 		void        GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
 		void        GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");};
@@ -138,5 +136,4 @@
 		void        GetInputValue(IssmDouble* pvalue,int enum_type){_error_("not implemented yet");};
 		void        GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type){_error_("not implemented yet");};
-		void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
 		Node*       GetNode(int node_number){_error_("Not implemented");};
 		IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
@@ -173,5 +170,5 @@
 		void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};
 		void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){_error_("not implemented yet");};
-		void        InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum){_error_("not implemented yet");};
+		void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
 		void        InputDuplicate(int original_enum,int new_enum){_error_("not implemented yet");};
 		void        InputScale(int enum_type,IssmDouble scale_factor){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17275)
@@ -178,11 +178,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::AddMaterialInput{{{*/
-void  Tria::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
-
-	_assert_(this->material);
-	this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
-}
-/*}}}*/
 /*FUNCTION Tria::CharacteristicLength{{{*/
 IssmDouble Tria::CharacteristicLength(void){
@@ -1018,9 +1011,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetMaterialInput(int inputenum) {{{*/
-Input* Tria::GetMaterialInput(int inputenum){
-	return this->material->inputs->GetInput(inputenum);
-}
-/*}}}*/
 /*FUNCTION Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype) {{{*/
 void Tria::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){
@@ -1173,17 +1161,4 @@
 
 }/*}}}*/
-/*FUNCTION Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
-void Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
-
-	Input* input=this->material->inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
-
-	GaussTria* gauss=new GaussTria();
-	gauss->GaussVertex(this->GetNodeIndex(node));
-
-	input->GetInputValue(pvalue,gauss);
-	delete gauss;
-}
-/*}}}*/
 Node* Tria::GetNode(int node_number){/*{{{*/
 	_assert_(node_number>=0); 
@@ -1247,5 +1222,5 @@
 /*}}}*/
 /*FUNCTION Tria::InputDepthAverageAtBase {{{*/
-void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
+void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){
 
 	/*New input*/
@@ -1254,10 +1229,5 @@
 
 	/*copy input of enum_type*/
-	if (object_enum==MeshElementsEnum)
-	 oldinput=(Input*)this->inputs->GetInput(enum_type);
-	else if (object_enum==MaterialsEnum)
-	 oldinput=(Input*)this->material->inputs->GetInput(enum_type);
-	else
-	 _error_("object " << EnumToStringx(object_enum) << " not supported yet");
+	oldinput=(Input*)this->inputs->GetInput(enum_type);
 	if(!oldinput)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
 	newinput=(Input*)oldinput->copy();
@@ -1267,10 +1237,5 @@
 
 	/*Add new input to current element*/
-	if (object_enum==MeshElementsEnum)
-	 this->inputs->AddInput((Input*)newinput);
-	else if (object_enum==MaterialsEnum)
-	 this->material->inputs->AddInput((Input*)newinput);
-	else
-	 _error_("object " << EnumToStringx(object_enum) << " not supported yet");
+	this->inputs->AddInput((Input*)newinput);
 }
 /*}}}*/
@@ -1408,7 +1373,20 @@
 					}
 					break;
-				/*Material will take care of it*/ 
-				case MaterialsRheologyBbarEnum:break;
-				case DamageDbarEnum: break;
+				case MaterialsRheologyBbarEnum:
+					if(iomodel->Data(MaterialsRheologyBEnum)){
+						for(j=0;j<3;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[tria_vertex_ids[j]-1];
+						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
+						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
+						this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+					}
+					break;
+				case DamageDbarEnum:
+					if(iomodel->Data(DamageDEnum)){
+						for(j=0;j<3;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[tria_vertex_ids[j]-1];
+						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
+						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
+						this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+					}
+					break;
 				default:
 					_error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
@@ -1476,10 +1454,5 @@
 		}
 		/*update input*/
-		if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==DamageDEnum || name==DamageDbarEnum){
-			material->inputs->AddInput(new TriaInput(name,values,P1Enum));
-		}
-		else{
-			this->inputs->AddInput(new TriaInput(name,values,P1Enum));
-		}
+		this->inputs->AddInput(new TriaInput(name,values,P1Enum));
 		break;
 
@@ -1490,10 +1463,5 @@
 		}
 		/*update input*/
-		if(name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==DamageDEnum || name==DamageDbarEnum){
-			material->inputs->AddInput(new TriaInput(name,values,P1Enum));
-		}
-		else{
-			this->inputs->AddInput(new TriaInput(name,values,P1Enum));
-		}
+		this->inputs->AddInput(new TriaInput(name,values,P1Enum));
 		break;
 
@@ -1603,61 +1571,4 @@
 	}
 
-}
-/*}}}*/
-/*FUNCTION Tria::IsInput{{{*/
-bool Tria::IsInput(int name){
-	if (
-				name==ThicknessEnum ||
-				name==SurfaceEnum ||
-				name==BathymetryEnum ||
-				name==BedEnum ||
-				name==MaskGroundediceLevelsetEnum ||
-				name==MaskIceLevelsetEnum ||
-				name==IceMaskNodeActivationEnum ||
-				name==LevelsetfunctionSlopeXEnum ||
-				name==LevelsetfunctionSlopeYEnum ||
-				name==SurfaceSlopeXEnum ||
-				name==SurfaceSlopeYEnum ||
-				name==BasalforcingsMeltingRateEnum ||
-				name==WatercolumnEnum || 
-				name==SurfaceforcingsMassBalanceEnum ||
-				name==SurfaceAreaEnum||
-				name==VxEnum ||
-				name==VyEnum ||
-				name==PressureEnum ||
-				name==InversionVxObsEnum ||
-				name==InversionVyObsEnum ||
-				name==FrictionCoefficientEnum ||
-				name==MaterialsRheologyBbarEnum ||
-				name==DamageDbarEnum ||
-				name==GradientEnum ||
-				name==OldGradientEnum ||
-				name==ConvergedEnum ||
-				name==SedimentHeadOldEnum ||
-				name==SedimentHeadEnum ||
-				name==EplHeadEnum ||
-				name==EplHeadOldEnum ||
-				name==HydrologydcEplThicknessOldEnum ||
-				name==HydrologydcEplInitialThicknessEnum ||
-				name==HydrologydcEplThicknessEnum ||
-				name==HydrologydcMaskEplactiveNodeEnum ||
-				name==MeshVertexonbedEnum ||
-				name==WaterTransferEnum ||
-				name==QmuVxEnum ||
-				name==QmuVyEnum ||
-				name==QmuPressureEnum ||
-				name==QmuBedEnum ||
-				name==QmuThicknessEnum ||
-				name==QmuSurfaceEnum ||
-				name==QmuTemperatureEnum ||
-				name==QmuMeltingEnum || 
-				name==QmuMaskGroundediceLevelsetEnum ||
-				name==QmuMaskIceLevelsetEnum ||
-				name==GiaWEnum || 
-				name==GiadWdtEnum
-		){
-		return true;
-	}
-	else return false;
 }
 /*}}}*/
@@ -2299,7 +2210,5 @@
 
 	/*Spawn material*/
-	seg->material=(Material*)this->material->copy();
-	delete seg->material->inputs;
-	seg->material->inputs=(Inputs*)this->material->inputs->SpawnSegInputs(index1,index2);
+	seg->material=(Material*)this->material->copy(seg);
 
 	/*recover nodes, material and matpar: */
@@ -3297,12 +3206,5 @@
 
 	for(int i=0;i<num_controls;i++){
-
-		if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==DamageDbarEnum){
-			input=(Input*)material->inputs->GetInput(control_type[i]); _assert_(input);
-		}
-		else{
-			input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
-		}
-
+		input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
 		if (input->ObjectEnum()!=ControlInputEnum){
 			_error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
@@ -3323,12 +3225,6 @@
 
 	int vertexpidlist[NUMVERTICES];
-	Input* input=NULL;
-
-	if(enum_type==MaterialsRheologyBbarEnum || enum_type==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(enum_type);
-	}
-	else{
-		input=inputs->GetInput(enum_type);
-	}
+
+	Input* input=inputs->GetInput(enum_type);
 	if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
 	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
@@ -3341,12 +3237,5 @@
 void Tria::ControlInputScaleGradient(int enum_type,IssmDouble scale){
 
-	Input* input=NULL;
-
-	if(enum_type==MaterialsRheologyBbarEnum || enum_type==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(enum_type);
-	}
-	else{
-		input=inputs->GetInput(enum_type);
-	}
+	Input* input=inputs->GetInput(enum_type);
 	if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
 	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
@@ -3360,12 +3249,6 @@
 	IssmDouble grad_list[NUMVERTICES];
 	Input* grad_input=NULL;
-	Input* input=NULL;
-
-	if(enum_type==MaterialsRheologyBbarEnum || enum_type==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(enum_type);
-	}
-	else{
-		input=inputs->GetInput(enum_type);
-	}
+
+	Input* input=inputs->GetInput(enum_type);
 	if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
 	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
@@ -3381,11 +3264,5 @@
 void Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){
 
-	Input* input=NULL;
-	if(control_enum==MaterialsRheologyBbarEnum || control_enum==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(control_enum);
-	}
-	else{
-		input=inputs->GetInput(control_enum);
-	}
+	Input* input=inputs->GetInput(control_enum);
 	if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
 	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
@@ -3499,5 +3376,5 @@
 	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	GradientIndexing(&vertexpidlist[0],control_index);
-	Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
+	Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
 	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
 
@@ -3548,5 +3425,5 @@
 	Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
 	Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
-	Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
+	Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3605,5 +3482,5 @@
 	Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
 	Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
-	Input* rheologyd_input=material->inputs->GetInput(DamageDbarEnum); _assert_(rheologyd_input);
+	Input* rheologyd_input=inputs->GetInput(DamageDbarEnum); _assert_(rheologyd_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4011,5 +3888,5 @@
 	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	Input* weights_input  =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);              _assert_(weights_input);
-	Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
+	Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
 
 	/* Start looping on the number of gaussian points: */
@@ -4573,10 +4450,5 @@
 
 	/*Get input (either in element or material)*/
-	if(control_enum==MaterialsRheologyBbarEnum || control_enum==DamageDbarEnum){
-		input=(Input*)this->material->inputs->GetInput(control_enum); _assert_(input);
-	}
-	else{
-		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
-	}
+	input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
 
 	/*Check that it is a ControlInput*/
@@ -4608,11 +4480,5 @@
 	new_input = new TriaInput(control_enum,values,P1Enum);
 
-	if(control_enum==MaterialsRheologyBbarEnum || control_enum==DamageDbarEnum){
-		input=(Input*)material->inputs->GetInput(control_enum); _assert_(input);
-	}
-	else{
-		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
-	}
-
+	input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
 	if (input->ObjectEnum()!=ControlInputEnum){
 		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
@@ -5150,4 +5016,7 @@
 
 					break;
+				case MaterialsRheologyBEnum:
+					this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values,P1Enum));
+					break;
 				default:
 					this->inputs->AddInput(new TriaInput(name,values,P1Enum));
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17275)
@@ -109,5 +109,5 @@
 		void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
 		void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
-		void        InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
+		void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
 		void        InputDuplicate(int original_enum,int new_enum);
 		void        InputScale(int enum_type,IssmDouble scale_factor);
@@ -202,5 +202,4 @@
 		void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
-		void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
 		IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
@@ -216,5 +215,4 @@
 		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
 		void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
-		Input*         GetMaterialInput(int inputenum);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
@@ -232,5 +230,4 @@
 		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
 		void	         InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
-		bool	         IsInput(int name);
 		void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
 		void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17275)
@@ -706,5 +706,5 @@
 
 	//First recover damage  using the element: */
-	element->GetMaterialInputValue(&damage,node,DamageDbarEnum);
+	element->GetInputValue(&damage,node,DamageDbarEnum);
 
 	//Recover our data:
Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 17275)
@@ -12,4 +12,5 @@
 #include "../../datastructures/datastructures.h"
 #include "../Update.h"
+class Element;
 class Elements;
 /*}}}*/
@@ -18,12 +19,10 @@
 
 	public: 
-		Inputs*  inputs;
 		virtual ~Material(){};
 		/*WARNING: input should not be public but it is an easy way to update B from T (using UpdateFromSolution) from Pentas*/
 
 		/*Numerics*/
-		virtual void       InputDuplicate(int original_enum,int new_enum)=0;
+		virtual Material*  copy(Element* element)=0;
 		virtual void       Configure(Elements* elements)=0;
-		virtual void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum)=0;
 		virtual void       GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
 		virtual void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
Index: /issm/trunk-jpl/src/c/classes/Materials/Materials.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Materials.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Materials.cpp	(revision 17275)
@@ -47,10 +47,2 @@
 }
 /*}}}*/
-/*FUNCTION Materials::InputDuplicate{{{*/
-void Materials::InputDuplicate(int input_enum,int output_enum){
-
-	for(int i=0;i<this->Size();i++){
-		Material* material=(Material*)this->GetObjectByOffset(i);
-		material->InputDuplicate(input_enum,output_enum);
-	}
-}
Index: /issm/trunk-jpl/src/c/classes/Materials/Materials.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Materials.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Materials.h	(revision 17275)
@@ -24,5 +24,4 @@
 		/*numerics*/
 		void  Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
-		void  InputDuplicate(int input_enum,int output_enum);
 
 };
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 17275)
@@ -29,6 +29,6 @@
 /*FUNCTION Matice::Matice(){{{*/
 Matice::Matice(){
-	this->inputs=NULL;
 	this->helement=NULL;
+	this->element=NULL;
 	return;
 }
@@ -42,14 +42,9 @@
 	/*Initialize id*/
 	this->mid=matice_mid;
-
-	/*Initialize inputs*/
-	this->inputs=new Inputs();
-
-	/*Initialize inputs from IoModel*/
-	this->InputUpdateFromIoModel(index,iomodel);
 
 	/*Hooks: */
 	matice_eid=index+1;
 	this->helement=new Hook(&matice_eid,1);
+	this->element=NULL;
 
 	return;
@@ -60,5 +55,4 @@
 Matice::~Matice(){
 	delete helement;
-	delete inputs;
 	return;
 }
@@ -71,6 +65,4 @@
 	_printf_("Matice:\n");
 	_printf_("   mid: " << mid << "\n");
-	_printf_("   inputs:\n");
-	inputs->Echo();
 	_printf_("   element:\n");
 	helement->Echo();
@@ -82,6 +74,4 @@
 	_printf_("Matice:\n");
 	_printf_("   mid: " << mid << "\n");
-	_printf_("   inputs:\n");
-	inputs->DeepEcho();
 	_printf_("   element:\n");
 	helement->Echo();
@@ -98,5 +88,5 @@
 }
 /*}}}*/
-/*FUNCTION Matice::copy {{{*/
+/*FUNCTION Matice::copy() {{{*/
 Object* Matice::copy() {
 
@@ -110,6 +100,22 @@
 	matice->mid=this->mid;
 	matice->helement=(Hook*)this->helement->copy();
-	if(this->inputs) matice->inputs=(Inputs*)this->inputs->Copy();
-	else  matice->inputs=new Inputs();
+	matice->element =(Element*)this->helement->delivers();
+
+	return matice;
+}
+/*}}}*/
+/*FUNCTION Matice::copy(Element* element) {{{*/
+Material* Matice::copy(Element* element_in) {
+
+	/*Output*/
+	Matice* matice=NULL;
+
+	/*Initialize output*/
+	matice=new Matice();
+
+	/*copy fields: */
+	matice->mid=this->mid;
+	matice->helement=(Hook*)this->helement->copy();
+	matice->element =element_in;
 
 	return matice;
@@ -124,4 +130,5 @@
 	 * datasets, using internal ids and offsets hidden in hooks: */
 	helement->configure((DataSet*)elementsin);
+	this->element  = (Element*)helement->delivers();
 }
 /*}}}*/
@@ -139,5 +146,5 @@
 	IssmDouble B,n;
 
-	inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
+	element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
 	n=this->GetN();
 
@@ -153,5 +160,5 @@
 	IssmDouble B,n;
 
-	inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
+	element->inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
 	n=this->GetN();
 
@@ -165,5 +172,5 @@
 	IssmDouble B;
 
-	inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
+	element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
 	return B;
 }
@@ -175,5 +182,5 @@
 	IssmDouble Bbar;
 
-	inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
+	element->inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
 	return Bbar;
 }
@@ -185,5 +192,5 @@
 	IssmDouble n;
 
-	inputs->GetInputAverage(&n,MaterialsRheologyNEnum);
+	element->inputs->GetInputAverage(&n,MaterialsRheologyNEnum);
 	return n;
 }
@@ -195,5 +202,5 @@
 	IssmDouble D;
 
-	inputs->GetInputAverage(&D,DamageDEnum);
+	element->inputs->GetInputAverage(&D,DamageDEnum);
 	return D;
 }
@@ -204,38 +211,6 @@
 	/*Output*/
 	IssmDouble Dbar;
-	inputs->GetInputAverage(&Dbar,DamageDbarEnum);
+	element->inputs->GetInputAverage(&Dbar,DamageDbarEnum);
 	return Dbar;
-}
-/*}}}*/
-/*FUNCTION Matice::GetVectorFromInputs{{{*/
-void  Matice::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){
-
-	/*Intermediaries*/
-	Element *element= NULL;
-
-	/*Recover element*/
-	element=(Element*)helement->delivers();
-
-	/*Check that input_enum is a material input*/
-	if (!IsInput(input_enum)) return;
-
-	switch(element->ObjectEnum()){
-
-		case TriaEnum:{
-
-							  /*Prepare index list*/
-							  int vertexpidlist[3];
-							  ((Tria*)element)->GetVertexPidList(&vertexpidlist[0]);
-
-								 /*Get input (either in element or material)*/
-								 Input* input=inputs->GetInput(input_enum);
-							  if(!input) _error_("Input " << EnumToStringx(input_enum) << " not found in material");
-
-							  /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
-							  input->GetVectorFromInputs(vector,&vertexpidlist[0]);}
-						break;
-
-		default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
-	}
 }
 /*}}}*/
@@ -546,100 +521,11 @@
 }
 /*}}}*/
-/*FUNCTION Matice::InputDuplicate{{{*/
-void  Matice::InputDuplicate(int original_enum,int new_enum){
-
-	/*Call inputs method*/
-	if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
-
-}
-/*}}}*/
 /*FUNCTION Matice::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
 void  Matice::InputUpdateFromVector(IssmDouble* vector, int name, int type){
 
-	/*Intermediaries*/
-	Element *element      = NULL;
-
-	/*Recover element*/
-	element=(Element*)helement->delivers();
-
-	/*Check that name is an element input*/
-	if (!IsInput(name)) return;
-
-	switch(type){
-
-		case VertexEnum:
-
-			switch(element->ObjectEnum()){
-
-				case TriaEnum: {
-					IssmDouble values[3];
-					for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()];
-					this->inputs->AddInput(new TriaInput(name,values,P1Enum));
-					return;
-				}
-				case PentaEnum: {
-					IssmDouble valuesp[6];
-					for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector
-					this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum));
-					return;
-				}
-				default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
-			}
-		default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
-	}
 }
 /*}}}*/
 /*FUNCTION Matice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
 void  Matice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
-
-	/*Intermediaries*/
-	Element    *element    = NULL;
-	Parameters *parameters = NULL;
-	int         meshtype;
-
-	/*Recover element*/
-	element=(Element*)helement->delivers();
-
-	/*Check that name is an element input*/
-	if (!IsInput(name)) return;
-
-	switch(type){
-
-		case VertexEnum:
-
-			switch(element->ObjectEnum()){
-
-				case TriaEnum: {
-					IssmDouble values[3];
-					for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //index into serial oriented vector 
-					this->inputs->AddInput(new TriaInput(name,values,P1Enum));
-					/*Special case for rheology B in 2D: Pourave land for this solution{{{*/
-					if(name==MaterialsRheologyBEnum){
-						/*Are we in 2D?:*/
-						if(element->ObjectEnum()==TriaEnum){
-							parameters=((Tria*)(element))->parameters;
-						}
-						else{
-							parameters=((Penta*)(element))->parameters;
-						}
-						parameters->FindParam(&meshtype,MeshTypeEnum);
-						if(meshtype==Mesh2DhorizontalEnum){
-							/*Duplicate rheology input: */
-							this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values,P1Enum));
-						}
-					}
-					/*}}}*/
-					return;
-				}
-				case PentaEnum:{
-					IssmDouble valuesp[6];
-					for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector
-					this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum));
-					return;
-				}
-				default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
-			}
-		default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
-	}
 
 }
@@ -665,167 +551,2 @@
 }
 /*}}}*/
-/*FUNCTION Matice::InputUpdateFromIoModel{{{*/
-void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){
-
-	int i,j;
-	bool   control_analysis;
-	int    num_control_type;
-
-	/*Fetch parameters: */
-	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
-	if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
-
-	if(iomodel->meshtype==Mesh2DhorizontalEnum){
-
-		/*Intermediaries*/
-		const int  num_vertices = 3; //Tria has 3 vertices
-		IssmDouble nodeinputs[num_vertices];
-		IssmDouble cmmininputs[num_vertices];
-		IssmDouble cmmaxinputs[num_vertices];
-
-		/*Get B*/
-		if (iomodel->Data(MaterialsRheologyBEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs,P1Enum));
-		}
-
-		/*Get n*/
-		if (iomodel->Data(MaterialsRheologyNEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
-			this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum));
-		}
-
-		/*Get D:*/
-		if (iomodel->Data(DamageDEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum));
-		}
-
-		/*Control Inputs*/
-		if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
-			for(i=0;i<num_control_type;i++){
-				switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
-					case MaterialsRheologyBbarEnum:
-						if (iomodel->Data(MaterialsRheologyBEnum)){
-							_assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
-							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+j]-1];
-							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
-						}
-						break;
-					case DamageDbarEnum:
-						if (iomodel->Data(DamageDEnum)){
-							_assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
-							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1];
-							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
-						}
-						break;
-				}
-			}
-		}
-	}
-	else if(iomodel->meshtype==Mesh2DverticalEnum){
-
-		/*Intermediaries*/
-		const int num_vertices = 3; //Tria has 3 vertices
-		IssmDouble    nodeinputs[num_vertices];
-		IssmDouble    cmmininputs[num_vertices];
-		IssmDouble    cmmaxinputs[num_vertices];
-
-		/*Get B*/
-		if (iomodel->Data(MaterialsRheologyBEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new TriaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum));
-		}
-
-		/*Get n*/
-		if (iomodel->Data(MaterialsRheologyNEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
-			this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum));
-		}
-
-		/*Get D:*/
-		if (iomodel->Data(DamageDEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new TriaInput(DamageDEnum,nodeinputs,P1Enum));
-		}
-	}
-	else if(iomodel->meshtype==Mesh3DEnum){
-
-		/*Intermediaries*/
-		const int  num_vertices = 6; //Penta has 6 vertices
-		IssmDouble nodeinputs[num_vertices];
-		IssmDouble cmmininputs[num_vertices];
-		IssmDouble cmmaxinputs[num_vertices];
-
-		/*Get B*/
-		if (iomodel->Data(MaterialsRheologyBEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum));
-		}
-
-		/*Get n*/
-		if (iomodel->Data(MaterialsRheologyNEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
-			this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum));
-		}
-		
-		/*Get D*/
-		if (iomodel->Data(DamageDEnum)) {
-			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new PentaInput(DamageDEnum,nodeinputs,P1Enum));
-		}
-		
-		/*Control Inputs*/
-		if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
-			for(i=0;i<num_control_type;i++){
-				switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
-					case MaterialsRheologyBbarEnum:
-						if (iomodel->Data(MaterialsRheologyBEnum)){
-							_assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
-							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[iomodel->elements[num_vertices*index+j]-1];
-							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
-						}
-						break;
-					case DamageDbarEnum:
-						if (iomodel->Data(DamageDEnum)){
-							_assert_(iomodel->Data(DamageDEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
-							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+j]-1];
-							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(iomodel->elements[num_vertices*index+j]-1)*num_control_type+i];
-							this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
-						}
-						break;
-
-				}
-			}
-		}
-	}
-	else{
-		_error_("Mesh type not supported yet!");
-	}
-
-	return;
-}
-/*}}}*/
-/*FUNCTION Matice::IsInput{{{*/
-bool Matice::IsInput(int name){
-	if (
-				name==MaterialsRheologyBEnum ||
-				name==MaterialsRheologyBbarEnum ||
-				name==MaterialsRheologyNEnum ||
-				name==QmuMaterialsRheologyBEnum ||
-				name==DamageDEnum ||
-				name==DamageDbarEnum
-
-		){
-		return true;
-	}
-	else return false;
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 17275)
@@ -12,4 +12,5 @@
 class IoModel;
 class Elements;
+class Element;
 class Loads;
 class Nodes;
@@ -22,6 +23,7 @@
 
 	private: 
-		int   mid;
-		Hook *helement;
+		int      mid;
+		Hook    *helement;
+		Element *element;
 
 	public:
@@ -45,10 +47,9 @@
 		void  InputUpdateFromConstant(int constant, int name);
 		void  InputUpdateFromConstant(bool constant, int name);
-		void  InputUpdateFromIoModel(int index, IoModel* iomodel);
+		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented");};
 		/*}}}*/
 		/*Material virtual functions resolution: {{{*/
-		void   InputDuplicate(int original_enum,int new_enum);
 		void   Configure(Elements* elements);
-		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum);
+		Material*  copy(Element* element);
 		void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
@@ -66,5 +67,4 @@
 		IssmDouble GetDbar();
 		IssmDouble GetN();
-		bool       IsInput(int name);
 		/*}}}*/
 };
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 17275)
@@ -91,6 +91,4 @@
 	iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
 	iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
-
-	this->inputs=NULL; /*not used here*/
 }
 /*}}}*/
@@ -234,12 +232,4 @@
 
 /*Matpar management: */
-/*FUNCTION Matpar::InputDuplicate{{{*/
-void  Matpar::InputDuplicate(int original_enum,int new_enum){
-
-	   /*Call inputs method*/
-	   if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
-
-}
-/*}}}*/
 /*FUNCTION Matpar::Configure {{{*/
 void  Matpar::Configure(Elements* elementsin){
@@ -520,5 +510,4 @@
 }
 /*}}}*/
-
 /*FUNCTION Matpar::GetLithosphereShearModulus {{{*/			 
 IssmDouble Matpar::GetLithosphereShearModulus(){		 
@@ -541,7 +530,2 @@
 }		 
 /*}}}*/ 
-/*FUNCTION Matpar::IsInput{{{*/
-bool Matpar::IsInput(int name){
-	return false;
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 17274)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 17275)
@@ -76,10 +76,9 @@
 		void   InputUpdateFromConstant(int constant, int name);
 		void   InputUpdateFromConstant(bool constant, int name);
-		void   InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void   InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented");};
 		/*}}}*/
 		/*Material virtual functions resolution: {{{*/
-		void       InputDuplicate(int original_enum,int new_enum);
+		Material*  copy(Element* element){_error_("not implemented");};
 		void       Configure(Elements* elements);
-		void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
 		void       GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
 		void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
@@ -135,5 +134,4 @@
 		IssmDouble GetS0p(); 
 		IssmDouble GetMaterialParameter(int in_enum); 
-		bool       IsInput(int name);
 		/*}}}*/
 
Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 17275)
@@ -21,10 +21,4 @@
 			Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
 			element->GetVectorFromInputs(vector,name);
-		}
-
-		/*Look up in materials*/
-		for(i=0;i<femmodel->materials->Size();i++){
-			Material* material=(Material*)femmodel->materials->GetObjectByOffset(i);
-			material->GetVectorFromInputs(vector,name);
 		}
 	}
Index: /issm/trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp	(revision 17275)
@@ -14,8 +14,3 @@
 		element->InputDuplicate(original_enum,new_enum);
 	}
-	for(int i=0;i<femmodel->materials->Size();i++){
-		Material* material=(Material*)femmodel->materials->GetObjectByOffset(i);
-		material->InputDuplicate(original_enum,new_enum);
-	}
-
 }
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 17275)
@@ -45,5 +45,5 @@
 	}
 
-	/*Update elements and materials: */
+	/*Update elements: */
 	counter=0;
 	for (i=0;i<iomodel->numberofelements;i++){
@@ -51,7 +51,4 @@
 			element=(Element*)elements->GetObjectByOffset(counter);
 			element->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements.
-
-			material=(Material*)materials->GetObjectByOffset(counter);
-			material->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements.
 			counter++;
 		}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 17274)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 17275)
@@ -28,5 +28,4 @@
 	switch(iomodel->meshtype){
 		case Mesh2DhorizontalEnum:
-			materials->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum);
 			for(i=0;i<iomodel->numberofelements;i++){
 				if(iomodel->my_elements[i]) elements->AddObject(new Tria(i+1,i,i,iomodel,nummodels));
@@ -51,17 +50,19 @@
 	switch(materials_type){
 		case MaticeEnum:
-			iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,DamageDEnum);
+			iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
+			iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
+			iomodel->FetchDataToInput(elements,DamageDEnum);
 			for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
-			if(dakota_analysis){
-				switch(iomodel->meshtype){
-					case Mesh2DhorizontalEnum:
-						materials->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum);
-						break;
-					case Mesh3DEnum:
-						materials->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); 
-						break;
-					default:
-						_error_("Mesh not supported yet");
-				}
+			switch(iomodel->meshtype){
+				case Mesh2DhorizontalEnum:
+					elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+					elements->InputDuplicate(DamageDEnum,DamageDbarEnum);
+					if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum);
+					break;
+				case Mesh3DEnum:
+					if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); 
+					break;
+				default:
+					_error_("Mesh not supported yet");
 			}
 			break;