Index: /issm/trunk/src/c/Container/Inputs.cpp
===================================================================
--- /issm/trunk/src/c/Container/Inputs.cpp	(revision 4966)
+++ /issm/trunk/src/c/Container/Inputs.cpp	(revision 4967)
@@ -426,2 +426,20 @@
 }
 /*}}}*/
+/*FUNCTION Inputs::AXPY{{{1*/
+void  Inputs::AXPY(int YEnum, double scalar, int XEnum){
+	   
+	Input* xinput=NULL;
+	Input* yinput=NULL;
+
+	/*Find x and y inputs: */
+	xinput=(Input*)this->GetInput(XEnum);
+	yinput=(Input*)this->GetInput(YEnum);
+
+	/*some checks: */
+	if(!xinput) ISSMERROR(" input %s could not be found!",EnumAsString(XEnum));
+	if(!yinput) ISSMERROR(" input %s could not be found!",EnumAsString(YEnum));
+
+	/*Apply AXPY: */
+	yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
Index: /issm/trunk/src/c/Container/Inputs.h
===================================================================
--- /issm/trunk/src/c/Container/Inputs.h	(revision 4966)
+++ /issm/trunk/src/c/Container/Inputs.h	(revision 4967)
@@ -34,4 +34,5 @@
 		Input*  GetInput(int enum_name);
 		Inputs* SpawnTriaInputs(int* indices);
+		void    AXPY(int YEnum, double scalar, int XEnum);
 		
 		void GetParameterAverage(double* pvalue, int enum_type);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 4966)
+++ /issm/trunk/src/c/Makefile.am	(revision 4967)
@@ -435,4 +435,6 @@
 					./modules/InputControlConstrainx/InputControlConstrainx.h\
 					./modules/InputControlConstrainx/InputControlConstrainx.cpp\
+					./modules/InputControlUpdatex/InputControlUpdatex.h\
+					./modules/InputControlUpdatex/InputControlUpdatex.cpp\
 					./modules/SurfaceAreax/SurfaceAreax.h\
 					./modules/SurfaceAreax/SurfaceAreax.cpp\
@@ -978,4 +980,6 @@
 					./modules/InputControlConstrainx/InputControlConstrainx.h\
 					./modules/InputControlConstrainx/InputControlConstrainx.cpp\
+					./modules/InputControlUpdatex/InputControlUpdatex.h\
+					./modules/InputControlUpdatex/InputControlUpdatex.cpp\
 					./modules/SurfaceAreax/SurfaceAreax.h\
 					./modules/SurfaceAreax/SurfaceAreax.cpp\
Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp	(revision 4966)
+++ /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp	(revision 4967)
@@ -15,10 +15,8 @@
 	int  dim;
 	int  numberofvertices;
-	bool extrude;
 	Vec  gradient = NULL;
 	
 	/*retrieve some parameters: */
 	parameters->FindParam(&dim,DimEnum);
-	parameters->FindParam(&extrude,ExtrudeParamEnum);
 	numberofvertices=vertices->NumberOfVertices();
 
@@ -36,7 +34,4 @@
 	VecAssemblyEnd(gradient);
 
-	/*Extrude if needed: */
-	if(dim==3 && extrude) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0);
-
 	/*Assign output pointers: */
 	*pgradient=gradient;
Index: /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.cpp
===================================================================
--- /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.cpp	(revision 4967)
+++ /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.cpp	(revision 4967)
@@ -0,0 +1,18 @@
+/*!\file InputControlUpdatex
+ * \brief: Y=Y+aX operation on inputs.
+ */
+
+#include "./InputControlUpdatex.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void InputControlUpdatex(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,double scalar,bool save_parameter){
+
+	/*Go through elemnets, and ask to carry out the operation on inputs: */
+	for(int i=0;i<elements->Size();i++){
+		Element* element=(Element*)elements->GetObjectByOffset(i);
+		element->InputControlUpdate(scalar,save_parameter);
+	}
+}
Index: /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.h
===================================================================
--- /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.h	(revision 4967)
+++ /issm/trunk/src/c/modules/InputControlUpdatex/InputControlUpdatex.h	(revision 4967)
@@ -0,0 +1,13 @@
+/*!\file:  InputControlUpdatex.h
+ * \brief header file for field extrusion
+ */ 
+
+#ifndef _INPUTCONTROLUPDATEX_H
+#define _INPUTCONTROLUPDATEX_H
+
+#include "../../Container/Container.h"
+
+/* local prototypes: */
+void InputControlUpdatex(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters, double scalar,bool save_parameter);
+
+#endif 
Index: /issm/trunk/src/c/modules/InputDuplicatex/InputDuplicatex.cpp
===================================================================
--- /issm/trunk/src/c/modules/InputDuplicatex/InputDuplicatex.cpp	(revision 4966)
+++ /issm/trunk/src/c/modules/InputDuplicatex/InputDuplicatex.cpp	(revision 4967)
@@ -12,5 +12,4 @@
 	
 	/*Go through elemnets, and ask to reinitialie the input: */
-
 	int      i;
 	for(i=0;i<elements->Size();i++){
@@ -18,4 +17,8 @@
 		element->InputDuplicate(original_enum,new_enum);
 	}
+	for(i=0;i<materials->Size();i++){
+		Material* material=(Material*)materials->GetObjectByOffset(i);
+		//material->InputDuplicate(original_enum,new_enum);
+	}
 
 }
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 4966)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 4967)
@@ -40,4 +40,7 @@
 	IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
 	IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+	if(iomodel->control_analysis && strcmp(iomodel->control_type,"rheology_B")==0){
+		IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+	}
 	
 	/*Create elements and materials: */
@@ -61,4 +64,5 @@
 	xfree((void**)&iomodel->rheology_B);
 	xfree((void**)&iomodel->rheology_n);
+	xfree((void**)&iomodel->control_parameter);
 
 	/*Add new constrant material property tgo materials, at the end: */
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 4966)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 4967)
@@ -47,5 +47,7 @@
 		IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
 		IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
-		IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+		if(strcmp(iomodel->control_type,"drag_coefficient")==0){
+			IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+		}
 	}
 
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/UpdateElementsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/UpdateElementsDiagnosticStokes.cpp	(revision 4966)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/UpdateElementsDiagnosticStokes.cpp	(revision 4967)
@@ -48,5 +48,7 @@
 		IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
 		IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
-		IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+		if(strcmp(iomodel->control_type,"drag_coefficient")==0){
+			IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+		}
 	}
 
Index: /issm/trunk/src/c/modules/modules.h
===================================================================
--- /issm/trunk/src/c/modules/modules.h	(revision 4966)
+++ /issm/trunk/src/c/modules/modules.h	(revision 4967)
@@ -28,4 +28,5 @@
 #include "./InputAXPYx/InputAXPYx.h"
 #include "./InputControlConstrainx/InputControlConstrainx.h"
+#include "./InputControlUpdatex/InputControlUpdatex.h"
 #include "./InputConvergencex/InputConvergencex.h"
 #include "./InputDuplicatex/InputDuplicatex.h"
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 4967)
@@ -72,4 +72,5 @@
 		virtual void   GetVectorFromInputs(Vec vector,int NameEnum)=0;
 		virtual void   InputAXPY(int YEnum, double scalar, int XEnum)=0;
+		virtual void   InputControlUpdate(double scalar,bool save_parameter)=0;
 		virtual void   InputControlConstrain(int control_type,double cm_min, double cm_max)=0;
 		virtual bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4966)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4967)
@@ -1026,4 +1026,60 @@
 }
 /*}}}*/
+/*FUNCTION Penta::InputControlUpdate{{{1*/
+void  Penta::InputControlUpdate(double scalar,bool save_parameter){
+
+	/*Intermediary*/
+	Input* input=NULL;
+	double cm_min,cm_max;
+	int    control_type;
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&cm_min,CmMinEnum);
+	this->parameters->FindParam(&cm_max,CmMaxEnum);
+	this->parameters->FindParam(&control_type,ControlTypeEnum);
+
+
+	/*Rheology*/
+	if(control_type==RheologyB2dEnum){
+
+		/*First, get revert to previous parameter value (kept in ControlParameter input)*/
+		matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBEnum);
+
+		/*Now, update using the gradient new = old + scalar * gradient*/
+		//matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum);
+		// For now: Gradient is in element (TO BE CHANGED) and parameter in matice
+		Input* input_B   =(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input_B);
+		Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
+		input_B->AXPY(input_Grad,scalar);
+
+		/*Constrain constrain input*/
+		input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input);
+		input->Constrain(cm_min,cm_max);
+
+		/*Finally: save input if requested*/
+		if (save_parameter) matice->inputs->DuplicateInput(RheologyBEnum,ControlParameterEnum);
+
+	}
+	else if(control_type==DragCoefficientEnum){
+
+		/*First, get revert to previous parameter value (kept in ControlParameter input)*/
+		this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
+
+		/*Now, update using the gradient new = old + scalar * gradient*/
+		this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
+
+		/*Constrain input*/
+		input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
+		input->Constrain(cm_min,cm_max);
+
+		/*Finally: save input if requested*/
+		if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
+
+	}
+	else{
+		ISSMERROR("control type %s not implemented yet",EnumAsString(control_type));
+	}
+}
+/*}}}*/
 /*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
 void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
@@ -1156,11 +1212,11 @@
 
 	/*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
-	for (i=0;i<this->inputs->Size();i++){
-		input=(Input*)this->inputs->GetObjectByOffset(i);
-		if (input->EnumType()==enum_type){
-			found=true; break;
-		}
-	}
-	if (!found) ISSMERROR("Input %s not found in penta->inputs",EnumAsString(enum_type));
+	if (enum_type==RheologyB2dEnum){
+		input=this->matice->inputs->GetInput(enum_type);
+	}
+	else{
+		input=this->inputs->GetInput(enum_type);
+	}
+	if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
 
 	/*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 
@@ -1748,8 +1804,4 @@
 		this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
 
-	}
-	if (iomodel->rheology_B) {
-		for(i=0;i<6;i++)nodeinputs[i]=iomodel->rheology_B[penta_vertex_ids[i]-1];
-		this->inputs->AddInput(new PentaVertexInput(RheologyBEnum,nodeinputs));
 	}
 	if (iomodel->control_parameter) {
@@ -5392,6 +5444,5 @@
 				name==VzEnum ||
 				name==TemperatureEnum ||
-				name==RheologyBEnum ||
-				name==RheologyNEnum ||
+				name==ControlParameterEnum ||
 				name==FitEnum ||
 				name==DragCoefficientEnum ||
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4967)
@@ -87,4 +87,5 @@
 		void   GradjDrag(Vec gradient);
 		void   InputAXPY(int YEnum, double scalar, int XEnum);
+		void   InputControlUpdate(double scalar,bool save_parameter);
 		void   InputControlConstrain(int control_type,double cm_min, double cm_max);
 		bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4966)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4967)
@@ -416,5 +416,10 @@
 
 			/*update input*/
-			this->inputs->AddInput(new TriaVertexInput(name,values));
+			if (name==RheologyB2dEnum || name==RheologyBEnum){
+				matice->inputs->AddInput(new TriaVertexInput(name,values));
+			}
+			else{
+				this->inputs->AddInput(new TriaVertexInput(name,values));
+			}
 			return;
 
@@ -628,5 +633,5 @@
 		}
 		else if (control_type==RheologyB2dEnum){
-			inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
+			matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
 			Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
 		}
@@ -965,5 +970,5 @@
 	adjointx_input=inputs->GetInput(AdjointxEnum);
 	adjointy_input=inputs->GetInput(AdjointyEnum);
-	rheologyb_input=inputs->GetInput(RheologyB2dEnum);
+	rheologyb_input=matice->inputs->GetInput(RheologyB2dEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1229,7 +1234,58 @@
 	yinput->AXPY(xinput,scalar);
 
-	/*Move input to Material if required (needed if control method) TO BE IMPROVED*/
-	if (YEnum==RheologyBEnum || YEnum==RheologyB2dEnum){
-		this->matice->inputs->AddInput((Input*)yinput->copy());
+}
+/*}}}*/
+/*FUNCTION Tria::InputControlUpdate{{{1*/
+void  Tria::InputControlUpdate(double scalar,bool save_parameter){
+
+	/*Intermediary*/
+	Input* input=NULL;
+	double cm_min,cm_max;
+	int    control_type;
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&cm_min,CmMinEnum);
+	this->parameters->FindParam(&cm_max,CmMaxEnum);
+	this->parameters->FindParam(&control_type,ControlTypeEnum);
+
+	/*Rheology*/
+	if(control_type==RheologyB2dEnum){
+
+		/*First, get revert to previous parameter value (kept in ControlParameter input)*/
+		matice->inputs->DuplicateInput(ControlParameterEnum,RheologyB2dEnum);
+
+		/*Now, update using the gradient new = old + scalar * gradient*/
+		//matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum);
+		// For now: Gradient is in element (TO BE CHANGED) and parameter in matice
+		Input* input_B   =(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input_B);
+		Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
+		input_B->AXPY(input_Grad,scalar);
+
+		/*Constrain constrain input*/
+		input=(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input);
+		input->Constrain(cm_min,cm_max);
+
+		/*Finally: save input if requested*/
+		if (save_parameter) matice->inputs->DuplicateInput(RheologyB2dEnum,ControlParameterEnum);
+
+	}
+	else if(control_type==DragCoefficientEnum){
+
+		/*First, get revert to previous parameter value (kept in ControlParameter input)*/
+		this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
+
+		/*Now, update using the gradient new = old + scalar * gradient*/
+		this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
+
+		/*Constrain input*/
+		input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
+		input->Constrain(cm_min,cm_max);
+
+		/*Finally: save input if requested*/
+		if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
+
+	}
+	else{
+		ISSMERROR("control type %s not implemented yet",EnumAsString(control_type));
 	}
 
@@ -1319,5 +1375,5 @@
 
 	/*Call inputs method*/
-	inputs->DuplicateInput(original_enum,new_enum);
+	if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
 
 }
@@ -1344,12 +1400,11 @@
 
 	/*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
-	for (i=0;i<this->inputs->Size();i++){
-		input=(Input*)this->inputs->GetObjectByOffset(i);
-		if (input->EnumType()==enum_type){
-			found=true;
-			break;
-		}
-	}
-	if (!found) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
+	if (enum_type==RheologyB2dEnum){
+		input=this->matice->inputs->GetInput(enum_type);
+	}
+	else{
+		input=this->inputs->GetInput(enum_type);
+	}
+	if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
 
 	/*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 
@@ -2122,8 +2177,4 @@
 		if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
 		this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
-	}
-	if (iomodel->rheology_B) {
-		for(i=0;i<3;i++)nodeinputs[i]=iomodel->rheology_B[tria_vertex_ids[i]-1];
-		this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs));
 	}
 	if (iomodel->control_parameter) {
@@ -6044,9 +6095,7 @@
 				name==MeltingRateEnum ||
 				name==AccumulationRateEnum ||
+				name==ControlParameterEnum ||
 				name==VxEnum ||
 				name==VyEnum ||
-				name==RheologyBEnum ||
-				name==RheologyB2dEnum ||
-				name==RheologyNEnum ||
 				name==FitEnum ||
 				name==DragCoefficientEnum ||
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4967)
@@ -83,4 +83,5 @@
 		void   GradjDrag(Vec gradient);
 		void   InputAXPY(int YEnum, double scalar, int XEnum);
+		void   InputControlUpdate(double scalar,bool save_parameter);
 		void   InputControlConstrain(int control_type,double cm_min, double cm_max);
 		bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
Index: /issm/trunk/src/c/objects/Materials/Material.h
===================================================================
--- /issm/trunk/src/c/objects/Materials/Material.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Materials/Material.h	(revision 4967)
@@ -18,4 +18,8 @@
 	public: 
 		virtual       ~Material(){};
+
+		/*Numerics*/
+		virtual void   InputDuplicate(int original_enum,int new_enum)=0;
+
 };
 #endif
Index: /issm/trunk/src/c/objects/Materials/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 4966)
+++ /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 4967)
@@ -55,4 +55,10 @@
 			this->inputs->AddInput(new TriaVertexInput(RheologyNEnum,nodeinputs));
 		}
+
+		/*Get control_parameter*/
+		if (iomodel->control_parameter) {
+			for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)];
+			this->inputs->AddInput(new TriaVertexInput(ControlParameterEnum,nodeinputs));
+		}
 	}
 
@@ -74,4 +80,10 @@
 			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->rheology_n[index];
 			this->inputs->AddInput(new PentaVertexInput(RheologyNEnum,nodeinputs));
+		}
+
+		/*Get control_parameter*/
+		if (iomodel->control_parameter) {
+			for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)];
+			this->inputs->AddInput(new PentaVertexInput(ControlParameterEnum,nodeinputs));
 		}
 	}
@@ -531,7 +543,41 @@
 }
 /*}}}*/
+/*FUNCTION Matice::InputDuplicate{{{1*/
+void  Matice::InputDuplicate(int original_enum,int new_enum){
+
+	/*Call inputs method*/
+	printf("%s\n",EnumAsString(original_enum));
+	if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
+
+}
+/*}}}*/
 /*FUNCTION Matice::InputUpdateFromVector(double* vector, int name, int type) {{{1*/
 void  Matice::InputUpdateFromVector(double* vector, int name, int type){
-	/*Nothing updated yet*/
+
+	/*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->Enum()){
+
+				case TriaEnum:
+					double values[3];
+					for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexDof()];
+					this->inputs->AddInput(new TriaVertexInput(name,values));
+					return;
+
+				default: ISSMERROR("element %s not implemented yet",EnumAsString(element->Enum()));
+			}
+		default: ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
+	}
 }
 /*}}}*/
@@ -566,2 +612,15 @@
 }
 /*}}}*/
+/*FUNCTION Matice::IsInput{{{1*/
+bool Matice::IsInput(int name){
+	if (
+				name==RheologyBEnum ||
+				name==RheologyB2dEnum ||
+				name==RheologyNEnum ||
+				name==ControlParameterEnum
+		){
+		return true;
+	}
+	else return false;
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Materials/Matice.h
===================================================================
--- /issm/trunk/src/c/objects/Materials/Matice.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Materials/Matice.h	(revision 4967)
@@ -50,4 +50,7 @@
 		void  InputUpdateFromSolution(double* solution);
 		/*}}}*/
+		/*Material virtual functions resolution: {{{1*/
+		void   InputDuplicate(int original_enum,int new_enum);
+		/*}}}*/
 		/*Matice Numerics: {{{1*/
 		void   Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin);
@@ -60,4 +63,5 @@
 		double GetB2d();
 		double GetN();
+		bool   IsInput(int name);
 		/*}}}*/
 };
Index: /issm/trunk/src/c/objects/Materials/Matpar.h
===================================================================
--- /issm/trunk/src/c/objects/Materials/Matpar.h	(revision 4966)
+++ /issm/trunk/src/c/objects/Materials/Matpar.h	(revision 4967)
@@ -54,4 +54,7 @@
 		void   InputUpdateFromSolution(double* solution);
 		/*}}}*/
+		/*Material virtual functions resolution: {{{1*/
+		void   InputDuplicate(int original_enum,int new_enum){ISSMERROR("not implemented yet");};
+		/*}}}*/
 		/*Numerics: {{{1*/
 		double GetG();
Index: /issm/trunk/src/c/solutions/control_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/control_core.cpp	(revision 4966)
+++ /issm/trunk/src/c/solutions/control_core.cpp	(revision 4967)
@@ -23,6 +23,4 @@
 	double  eps_cm;
 	double  tolx;
-	double  cm_min;
-	double  cm_max;
 	bool    cm_gradient;
 	int     dim;
@@ -54,6 +52,4 @@
 	femmodel->parameters->FindParam(&eps_cm,EpsCmEnum);
 	femmodel->parameters->FindParam(&tolx,TolXEnum);
-	femmodel->parameters->FindParam(&cm_min,CmMinEnum);
-	femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
 	femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
 	femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
@@ -90,15 +86,6 @@
 		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
 
-		_printf_("%s\n","      updating parameter using optimized search scalar...");
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
-		InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum);
-
-		_printf_("%s\n","      constraining the new distribution...");    
-		InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
-		
-		_printf_("%s\n","      save new parameter...");
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum);
-		
-		_printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
+		_printf_("%s\n","      updating parameter using optimized search scalar..."); //true means update parameter and copy it onto ControlParameter input
+		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true);
 		
 		if(controlconvergence(J,fit,eps_cm,n)) break;
Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4966)
+++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4967)
@@ -32,6 +32,4 @@
 	double* optscal=NULL;
 	double* fit=NULL;
-	double  cm_min;
-	double  cm_max;
 	int     control_type;
 	int     analysis_type;
@@ -48,6 +46,4 @@
 	femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
 	femmodel->parameters->FindParam(&fit,NULL,FitEnum);
-	femmodel->parameters->FindParam(&cm_min,CmMinEnum);
-	femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
 	femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
 	femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
@@ -59,12 +55,6 @@
 	else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
 
-	/*Use ControlParameterEnum input to  reinitialize our input parameter: */
-	InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
-	
-	/*Use search scalar to shoot parameter in the gradient direction: */
-	InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum);
-
-	/*Constrain:*/
-	InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
+	/*update parameter according to scalar: */ //false means: do not copy updated parameter onto ControlParameter input
+	InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],false);
 
 	/*Run diagnostic with updated inputs: */
Index: /issm/trunk/src/m/solutions/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/control_core.m	(revision 4966)
+++ /issm/trunk/src/m/solutions/control_core.m	(revision 4967)
@@ -19,6 +19,4 @@
 	eps_cm=femmodel.parameters.EpsCm;
 	tolx=femmodel.parameters.TolX;
-	cm_min=femmodel.parameters.CmMin;
-	cm_max=femmodel.parameters.CmMax;
 	cm_gradient=femmodel.parameters.CmGradient;
 	control_steady=femmodel.parameters.ControlSteady;
@@ -55,13 +53,5 @@
 
 		displaystring('\n%s',['      updating parameter using optimized search scalar:']);
-		scalar=search_scalar*optscal(n);
-		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
-		[femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum);
-
-		displaystring('\n%s',['      constraning the new distribution...']);
-		[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
-
-		displaystring('\n%s',['      save new parameter...']);
-		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,ControlParameterEnum);
+		[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal(n),1);
 
 		disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);
Index: /issm/trunk/src/m/solutions/objectivefunctionC.m
===================================================================
--- /issm/trunk/src/m/solutions/objectivefunctionC.m	(revision 4966)
+++ /issm/trunk/src/m/solutions/objectivefunctionC.m	(revision 4967)
@@ -10,6 +10,4 @@
 analysis_type=femmodel.parameters.AnalysisType;
 isstokes=femmodel.parameters.IsStokes;
-cm_min=femmodel.parameters.CmMin;
-cm_max=femmodel.parameters.CmMax;
 
 %set current configuration
@@ -20,13 +18,6 @@
 end
 
-%Use ControlParameterEnum input to  reinitialize our input parameter:
-femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
-
 %Use search scalar to shoot parameter in the gradient direction:
-scalar=search_scalar*optscal;
-[femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum);
-
-%Constrain:
-[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
+[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal,0);
 
 %Run diagnostic with updated inputs:
Index: /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.cpp
===================================================================
--- /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.cpp	(revision 4967)
+++ /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.cpp	(revision 4967)
@@ -0,0 +1,68 @@
+/*\file InputControlUpdate.c
+*\brief: update elements properties using an input  vector
+*/
+
+#include "./InputControlUpdate.h"
+
+void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+/*input datasets: */
+Elements   *elements   = NULL;
+Nodes      *nodes      = NULL;
+Vertices   *vertices   = NULL;
+Loads      *loads      = NULL;
+Materials  *materials  = NULL;
+Parameters *parameters = NULL;
+double      scalar;
+double      update;
+
+/*Boot module: */
+MODULEBOOT();
+
+/*checks on arguments on the matlab side: */
+CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InputControlUpdateUsage);
+
+/*Input datasets: */
+FetchData((DataSet**)&elements,ELEMENTSIN);
+FetchData((DataSet**)&nodes,NODESIN);
+FetchData((DataSet**)&vertices,VERTICESIN);
+FetchData((DataSet**)&loads,LOADSIN);
+FetchData((DataSet**)&materials,MATERIALSIN);
+FetchParams(&parameters,PARAMETERSIN);
+FetchData(&scalar,SCALAR);
+FetchData(&update,UPDATE);
+
+/*configure: */
+elements->  Configure(elements,loads, nodes,vertices, materials,parameters);
+nodes->     Configure(elements,loads, nodes,vertices, materials,parameters);
+loads->     Configure(elements, loads, nodes,vertices, materials,parameters);
+
+/*call "x" code layer*/
+InputControlUpdatex(elements,nodes,vertices,loads, materials,parameters,scalar,(bool)update);
+
+/*write output datasets: */
+WriteData(ELEMENTS,elements);
+WriteData(NODES,nodes);
+WriteData(VERTICES,vertices);
+WriteData(LOADS,loads);
+WriteData(MATERIALS,materials);
+WriteParams(PARAMETERS,parameters);
+
+/*Free ressources: */
+delete elements;
+delete nodes;
+delete vertices;
+delete loads;
+delete materials;
+delete parameters;
+
+/*end module: */
+MODULEEND();
+}
+
+void InputControlUpdateUsage(void)
+{
+	_printf_("\n");
+	_printf_("   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,parameters,scalar);\n",__FUNCT__);
+	_printf_("\n");
+}
Index: /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.h
===================================================================
--- /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.h	(revision 4967)
+++ /issm/trunk/src/mex/InputControlUpdate/InputControlUpdate.h	(revision 4967)
@@ -0,0 +1,44 @@
+/*
+	InputControlUpdate.h
+*/
+
+#ifndef _INPUTCONTROLUPDATE_H
+#define _INPUTCONTROLUPDATE_H
+
+/* local prototypes: */
+void InputControlUpdateUsage(void);
+
+#include "../../c/modules/modules.h"
+#include "../../c/Container/Container.h"
+#include "../../c/shared/shared.h"
+#include "../../c/EnumDefinitions/EnumDefinitions.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "InputControlUpdate"
+
+/* serial input macros: */
+#define ELEMENTSIN (mxArray*)prhs[0]
+#define NODESIN (mxArray*)prhs[1]
+#define VERTICESIN (mxArray*)prhs[2]
+#define LOADSIN (mxArray*)prhs[3]
+#define MATERIALSIN (mxArray*)prhs[4]
+#define PARAMETERSIN (mxArray*)prhs[5]
+#define SCALAR (mxArray*)prhs[6]
+#define UPDATE (mxArray*)prhs[7]
+
+/* serial output macros: */
+#define ELEMENTS (mxArray**)&plhs[0]
+#define NODES (mxArray**)&plhs[1]
+#define VERTICES (mxArray**)&plhs[2]
+#define LOADS (mxArray**)&plhs[3]
+#define MATERIALS (mxArray**)&plhs[4]
+#define PARAMETERS (mxArray**)&plhs[5]
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  6
+#undef NRHS
+#define NRHS  8
+
+#endif  /* _INPUTCONTROLCONSTRAIN_H */
+
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 4966)
+++ /issm/trunk/src/mex/Makefile.am	(revision 4967)
@@ -22,4 +22,5 @@
 				ElementConnectivity\
 				InputControlConstrain \
+				InputControlUpdate \
 				InputConvergence\
 				GetSolutionFromInputs\
@@ -134,4 +135,7 @@
 			  InputControlConstrain/InputControlConstrain.h
 
+InputControlUpdate_SOURCES = InputControlUpdate/InputControlUpdate.cpp\
+										  InputControlUpdate/InputControlUpdate.h
+
 InputConvergence_SOURCES = InputConvergence/InputConvergence.cpp\
 			  InputConvergence/InputConvergence.h